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Section  1 


To  adequately  describe  the  strain-rate  sensitivity  and  creep  behavior  of 
P21S  during  isothermal  and  nonisothermal  load  conditions  requires  a  unified 
inelastic-strain  theory.  The  unified  inelastic-strain  theory  consisting  of  a  Bodner- 
Partom  flow  rule  [1]  with  modifications  can  adequately  model  the  response  of 
P21S.  Due  to  the  advanced  nature  of  these  theories,  many  finite  element 
progranas  do  not  include  these  formulations  directly  into  their  codes.  However, 
general  purpose  finite  element  programs,  such  as  ADINA6.0  [2},  do  allow  user- 
defined  subroutines  that  can  be  written  to  solve  unified  inelastic-strain  theories. 
This  manual  describes  this  unified  inelastic  strain  theory  for  (5215  and  provides 
the  associated  user-defined  subroutines  for  use  with  ADINA6.0.  The  manual  is 
designed  to  assist  the  user  with  implementing  these  subroutines  into  AD1NA6.0 
or  modifying  the  algorithms  within  these  subroutines  for  use  with  other  material 
behavior  models  or  finite  element  codes.  These  subroutines  and  associated  data 
files  given  in  the  appendices  are  available  by  electric  mail  or  on  a  magnetic 
medium  from  the  authors. 
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Section  2 

Bodner-Partom  Constitutive  Model 

This  section  presents  the  Bodner-Partom  constitutive  equations  and  the 
material  parameters  that  describe  the  P21S  response.  Several  formulations  of  the 
Bodner-Partom  model  can  characterize  inelastic  deformation  under  a  variety  of 
conditions  which  may  account  for  anisotropic,  isothermal,  and/or  non- 
isothermal  material  response.  The  formulation  presented  in  the  Section  2.1 
parallels  the  previous  isotropic  nonisothermal  theory  of  Chan,  Bodner,  and 
Lindholm  [3,4,5].  The  terminology  is  similar  to  that  of  Chan  and  Lindholm  [6]. 
Section  2.2  contains  the  material  parameters  that  characterize  the  P21S  material. 

2.1  Theory 


The  first  assumption  in  this  analysis  is  the  decomposition  of  the  total 
strain,  into  elastic,  thermal,  and  inelastic  components.  This  decomposition 
is  expressed  as: 


Eq.  1 


The  elastic  strains,  depend  on  the  current  stress  state,  elastic  modulus  E,  and 
Poisson's  ratio  v.  The  thermal  strain  components,  e||',  equal  the  product  of  the 
coefficient  of  thermal  expansion  and  the  difference  between  the  current  and 
reference  temperatures.  The  Bodner-Partom  flow  rule  governs  the  evolution  of 
the  inelastic  strains,  £-|^. 
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As  opposed  to  other  unified  inelastic  strain  formulations,  this  theory 
describes  the  directional  (kinematic)  hardening  with  a  special  directional 
hardening  term.  Other  theories  represent  directional  hardening  phenomena  with 
a  "back-stress"  modified  equivalent  stress  [7,8,9].  Introduction  of  the  directional 
hardening  term  alters  the  Bodner-Partom  flow  rule  by  replacing  the  previously 
known  variable  "drag-stress"  with  the  sum  of  isotropic  and  directional  hardening 
terms,  and  tP,  respectively.  These  two  hardening  terms  enter  into  the 
inelastic  strain  rate  equation  or  flow  law: 


Eq.2 


where  Dq  is  the  limiting  strain  rate,  sij  are  the  components  of  deviatoric  stress, 
and  J2  =  Sij  sy  /2. 


The  evolutions  of  and  Z^  have  similar  empirical  forms.  Each  equation 
consists  of  a  hardening  term,  a  thermal  recovery  term,  and  a  temperature  rate 
term.  The  isotropic  hardening  evolution  equation  with  these  three  terms  is: 

Z^  =  mi(Zi-zI)w‘*'-AiZi  —7-^ 

Eq.  3 


Zl-Z2 

azi^ 

f  Zi-zO 

^Z2 

ilZl-Z2j 

8T 

,Zi-Z2^ 

ar  J 

where  the  inelastic  work  rate  is  given  by: 


W 


Oij£!{. 


Eq.4 
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The  initial  isotopic  hardening,  Z^(0),  is  Zq.  The  material  parameters  associated 
with  the  isotopic  hardening  evolution  are  mi,  Zo,  Zi,  Z2,  Ai,  and  ri-  The 
thermal  differential  terms and  appropriately  scale  the  isotropic 

hardening  variable  when  inelastic  deformation  and  thermal  recovery  do  not 
occur  during  nonisothermal  conditions,  thus  preventing  Z^  from  passing  through 
maximum  or  minimum  values  with  temperature  changes.  The  treatment  of  these 
thermal  differential  terms  is  more  consistent  with  the  work  of  McDowell  [10]  and 
others  [11,12,13,14]  than  those  proposed  by  Chan,  Lindholm  and  Bodner  [4].  The 
thermal  differential  terms  appearing  in  the  user-defined  subroutines  of 
Appendices  A  and  B  reflect  the  fact  that  Zi  is  temperature  independent  for  P21S; 

azi 


J-e, 


3T 


=0. 


The  scalar  product  of  a  state  variable,  Pjj,  and  a  unit  stress  vector  uij 
yields  the  magnitude  of  the  directional  hardening  term: 


zD  =  PijUij, 


where: 


4ij  =  m2(Z3Uij-)3ij)wi^-A2Z^ 


Eq.5 


^^ijaZ3+  c  , 
vi|  +  Z — 


Z3 


dT 


and; 


vij  = 


Okl^kl 


^PklPkl 


Eq.7 


Eq.8 
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The  initial  directional  hardening  variable,  Z^^O)  is  set  to  zero.  The  material 

constants  associated  with  the  directional  hardening  evolution  equation  are  m2, 

dZ'\ 

Z3,  A2,  and  r2-  The  temperature  differential  term  — ~  appropriately  scales  the 

9T 

directional  hardening  variable  when  inelastic  deformation  and  thermal  recovery 
do  not  occur  during  nonisothermal  conditions.  In  particular,  without  these 
differential  terms  the  directional  hardening  accrued  at  one  temperature  may 
exceed  the  limiting  value  Z3  at  another  temperature,  which  is  not  physically 
realistic. 


Table  1  summarizes  the  15  material  parameters  that  characterize  the 
strain-rate  sensitivity  and  time-dependent  behavior  of  P21S.  The  number  of 
material  parameters  is  effectively  less  than  those  in  Table  1  after  applying  the 
usual  assumptions:  Ai  =  A2,  ri  =  r2,  and  Zq  =  Z2.  The  temperature  dependency 
of  these  constants  is  specific  to  |321S.  Other  materials,  such  as  Mar-M47  or 
B1900+Hf  [5],  require  different  temperature-dependent  parameters  to  represent 
their  particular  response. 
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Table  1 


Material  Parameters  in  Bodner-Partom  Constitutive  Model 


Parameter 

Units 

Description 

Temperature 

Dependency 

E 

MPa 

Elastic  modulus 

yes 

V 

— 

Poisson's  ratio 

no 

a 

1/°C 

Coefficient  of  thermal  expansion 

yes 

Do 

1 

C" 

w 

Limiting  shear  strain  rate 

no 

Zo 

MPa 

Initial  value  of  the  isotropic 
hardening  variable 

yes 

Zi 

MPa 

Limiting  (maximum)  value  of  Z 

no 

Z2 

MPa 

Fully  recovered  (minimum)  value 
ofZl 

yes 

Z3 

MPa 

Limiting  (maximum)  value  of 

yes 

mi 

MPa-1 

Hardening  rate  coefficient  of  Zl 

yes 

m2 

MPa-1 

Hardening  rate  coefficient  of  Z^ 

yes 

n 

— 

yes 

Ai 

s-i 

Recovery  coefficient  for  Z* 

yes 

A2 

s-i 

Recovery  coefficient  for  Zl^ 

yes 

ri 

— 

Recovery  exponent  for  Zl 

yes 

ri 

— 

Recovery  exponent  for  Zl^ 

yes 

2.2  B21S  Material  Parameters 

The  material  parameters  for  the  Bodner-Partom  model  with  directional 
hardening  were  determined  from  p21S  using  monotonic,  cyclic,  and  creep  test 
data  [15].  The  parameters  are  valid  for  a  wide  window  of  strain  rates  (10*3  to  10*^ 
1/s)  and  temperatures  (23°  -  815°C).  The  strategy  for  determining  the 
parameters  involves  a  number  of  steps.  First,  the  temperature-dependent 
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parameters  are  identified  and  the  values  of  the  temperature-independent 
parameters  are  estimated.  Then  at  each  temperature  where  experimental  data 
are  available,  the  temperature-dependent  parameters  are  determined  through  an 
iterative  process  to  minimize  the  differences  between  the  model  simulations  and 
experiments.  Mathematica  [16]  is  used  to  generate  the  model  simulations. 
Similar  to  other  inelastic  strain  theories,  the  set  of  material  parameter  for  any 
particular  load  case  is  not  unique.  Thus,  for  a  given  set  of  experimental  load 
responses,  a  range  of  values  are  suitable  for  each  material  parameter. 

The  resulting  set  of  temperature-dependent  parameters  becomes 
continuous  with  temperature  as  the  range  of  possible  values  is  decreased  for  each 
parameter.  The  response  can  be  very  sensitive  to  small  changes  in  some  of  the 
material  parameters  with  temperature,  especially  in  the  transition  regimes 
between  different  deformation  mechanisms.  For  p2IS,  a  transition  region  for 
inelastic  behavior  occurs  between  482°C  for  plasticity  and  at  650°C  for  power  law 
creep.  Thus,  an  anomalous  change  in  the  saturated  stress  level  can  occur  if  linear 
interpolation  of  material  parameters  is  used  within  this  transition  region. 

Since  no  experimental  data  are  available  within  this  transition  region, 
values  for  temperature-dependent  parameters  are  chosen  between  482°C  and 
650°C  so  that  the  resulting  saturated  stress  is  smooth  and  continuous  with 
temperature,  thus  reducing  the  anomalous  effects  of  linear  interpolation.  The 
final  version  of  the  material  parameters  (version  4.0)  for  (5215  are  shown  in  Table 
2. 
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Table  2 


Bodner-Partom  Material  Parameters  for  0215 

i 

i 


Temperatvire-Independent  Constants 


mi  =  0  MPa-^  (no  isotropic  hardening) 

Zi  =  1600  MPa 

ri  =  r2  =  3 

Do  =  1  X  104  s-l 

n  =  0.34 

Temperature-Dependent  Constants 


Temp. 

E 

a* 

n 

Zo  -  Z2 

Z3 

m2 

Ai  =  A2 

MPa 

10-6/‘’C 

MPa 

MPa 

MPa-1 

s-1 

23 

112000 

6.31 

4.8 

1550 

100 

0.35 

IHIIIIIQIIIIIII 

108000 

7.26 

3.5 

1300 

300 

035 

315 

4 

0 

0 

0 

390 

0 

0 

365 

0 

0 

0 

0 

500 

0 

0 

415 

4 

0 

0 

0 

660 

0 

0 

465 

<!• 

0 

0 

0 

960 

0 

0 

482 

98100 

8.15 

1,7 

1100 

1100 

5 

0.0076 

0 

1.5 

0 

1300 

0 

0 

525 

<t» 

0 

1.28 

0 

1670 

0 

0 

« 

0 

1.1 

0 

2100 

0 

0 

575 

4> 

0 

0.97 

0 

2600 

0 

0 

600 

0 

0 

0.82 

0 

3700 

10 

0 

650 

86600 

8.83 

0.74 

1000 

3800 

10 

0.21 

760 

77200 

9.27 

0.58 

600 

4000 

15 

1.0 

815 

72000 

9.49 

0.55 

300 

4100 

30 

2.0 

•  Secant  a  with  Tq  =  25°C 

0 '  Linear  interpolate  between  values  given  in  table. 

^  -  Use  functions  to  determine  values. 

Bold  were  values  determined  based  on  experiments. 

Italics  are  values  that  describe  smooth  and  continuous  saturated  stress  change. 


Functions 

E  =  112000  + 0.591  T- 0.061  T2  Tin°C 
a  =  6.22  X 10-6  +  4.01  x  lO'^  T  T  in  °C 
-  /-1.37xl04>, 

Ai  =  A2  =  5.8 X  10^  exp(  j  T in  °C 
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Section  3 


Installation  of  the  Subroutines 

Included  with  the  delivered  version  of  ADINA  6.0  are  overlay  files  that 
allow  users  to  incorporate  new  segments  of  code  into  certain  regions  within 
ADINA's  algorithms.  Two  particular  overlay  files,  ovl030u.f  and  ovl040u.f, 
contain  examples  of  user-defined  subroutines  that  incorporate  alternative 
material  response  models  into  ADINA.  Within  these  two  files,  the  user-defined 
subroutines,  CUSER2  and  CUSER3,  contain  sections  of  code  that  are  modified  to 
implement  the  Bodner-Partom  model  into  ADINA's  two-  and  three-dimensional 
elements,  respectively.  The  Bodner-Partom  versions  of  the  user-defined 
subroutines  for  the  two-dimensional  axisymmetric  and  three-dimensional  brick 
elements  reside  in  deliverable  files,  cuser2_dbeta.f  and  cuser3_dbeta.f, 
respectively  (see  Appendices  A  and  B).  The  last  step  of  the  installation  process  is 
to  integrate  the  subroutines  within  these  two  files  into  ADINA's  numerical 
algorithms. 

Integrating  the  Bodner-Partom  user-defined  subroutines  into  ADINA6.0 
can  be  accomplished  by  several  methods.  One  method  is  replacing  the 
subroutines  within  the  overlay  files  from  ADINA  with  those  of  cuser2_dbeta.f 
and  cuser3_dbeta.f.  Another  method  consists  of  linking  new  object  code  over  the 
pre-existing  object  code.  In  particular,  subroutines  of  files  cuser2_dbeta.f  and 
cuser3_dbeta.f  are  compiled  for  produce  two  object  files.  These  object  files  are 
then  linked  with  ADINA's  object  code  to  generate  a  new  executable  version  of 
ADINA  that  contains  the  Bodner-Partom  algorithms.  Unfortunately,  each  linker 
(or  loader)  is  machine  dependent.  For  example,  some  linkers  may  declare 
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multiple  subroutine  declarations  and  cease  the  operation.  Renaming  or  replacing 
the  subroutines  CUSER2  and  CUSER3  in  overlay  files,  ovl030u.f  and  ovl040u.f, 
respectively,  will  correct  these  errors.  Refer  to  your  FORTRAN??  manual  to 
determine  the  details  of  your  particular  linker  or  segment  loader. 

To  avoid  a  common  error,  check  the  precision  (single  or  double)  of  the 
existing  user-subroutines  found  in  ADDMA's  overlay  files  (e.g.,  olv030u.f).  Then, 
modify  the  variable  precision  of  the  subroutines  in  files  cuser2_dbeta.f  and 
cuser3_dbeta.f  to  be  of  the  same  precision  as  found  in  the  overlay  files. 
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Section  4 


Structure  and  Algorithms 

This  section  describes  the  structure  and  numerical  algorithms  of  the 
Bodner-Partom  modified  subroutines  which  characterize  the  strain-rate 
sensitivity  and  creep  behavior  of  P21S.  Several  tables  briefly  describe  subroutine 
nomenclatiu-e  and  variable  usage.  Then,  the  numerical  algorithm  that  solve  the 
Bodner-Partom  equations  is  outlined.  Hard  copies  of  the  subroutines  described 
in  this  section  reside  in  the  files  cuser2_dbeta.f  and  cuser3_dbeta.f  for  the  two- 
dimensional  axisymmetric  and  three-dimensional  brick  elements,  respectively, 
(see  Appendices  A  and  B).  Subroutine  structure  and  variables  usage  are  identical 
for  the  two-dimensional  and  three-dimensional  elements.  Thus  the  details  are 
only  presented  for  the  two-dimensional  case.  The  names  appended  with  2  and  3 
designate  the  two-  and  three-  dimensional  case,  respectively.  For  example, 
CUSER2  and  CUSER3  are  user-defined  subroutines  for  the  two-  and  three- 
dimensional  cases,  respectively. 

4.1  Subroutines  and  Common  Blocks 

The  functionality  of  subroutines  found  in  cuser2_dbeta.f  are  summarized 
in  Table  3.  Subroutines  lUSER  and  USER2  appear  in  the  overlay  file  ovl030u.f 
which  is  supplied  with  the  ADINA6.0  [2].  Subroutine  CUSER2  has  four  major 
functions  that  integrate  new  material  behavior  models  into  the  ADINA  code. 
The  four  functions  include  -- 1)  the  initialization  of  internal  state  variables,  2)  the 
determination  of  new  states  from  previous  strains  and  current  incremental 
strains,  3)  the  computation  of  incremental  stiffness  matrix,  D,  and  4)  output  of 


the  converged  solution.  A  minor  function  of  CUSER2  is  to  retrieve  and  store 
state  variables  from  ADINA's  working  variable,  ARRAY.  Subroutine  STGET2 
retrieves  while  STPUT2  stores  these  variables  from  ARRAY.  DBODNER2 
contain  the  numerical  algorithms  that  solve  the  Bodner-Partom  constitutive 
equations  found  in  Section  2.1 


Tables 


Subroutine  Names  and  Functions 


Subroutine* 

Called  From 

Function 

CUSER2 

IUSER2  and 
USER2 

to  determine  stresses  from  user- 
defined  material  behavior  models 

STGET2 

CUSER2 

to  retrieve  state  variables  from 
working  variable  ARRAY 

STPIJT2 

CUSER2 

to  store  state  variables  from  working 
variable  ARRAY 

DBODNER2 

CUSER2 

to  solve  the  constitutive  equations  of 
the  Bodner  Partom  model 

•Two-dimensional  subroutines  arc  listed.  The  three-dimensional  subroutines  are  designated  with 
a  suffix  3  rather  than  a  2. 


The  two  common  blocks  MATCONST  and  BPSTAT2  retain  values  of 
certain  variables  entering  from  other  subroutines.  MATCONST  transfers 
material  constants  and  associated  temperature  differentials  between  CUSER2 
and  DBODNER2.  Table  4  lists  variables  that  are  in  common  block  MATCONST. 
MATCONST  is  the  same  common  block  in  the  two-dimensional  element  as 
found  in  the  three-dimensional  elements.  Common  block  BPSTAT2  transfers  the 
values  of  the  internal  state  variables  (e.g.,  inelastic  strain,  isotropic  hardening, 
etc.)  between  all  the  subroutines.  Table  5  contains  the  state  variables  of  common 
block  BPSTAT2  and  BPSTAT3. 
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Table  4 


Internal  Variables  of  Common  Block  MATCONST 


Variable  Name 

Variable  Description 

E 

E 

elastic  modulus,  CTD(l) 

AN 

n 

kinetic  parameter,  CTD(2) 

zo;z2 

Z0,Z2 

fully  recovered  (miitimum)  value  of  isotropic 
hardening,  CTD(3) 

Z3 

Z3 

limiting  (maximum)  value  of  directional 
hardening,  CTD(4) 

AM2 

m2 

hardening  rate  coefficient  of  directional 
hardening,  CTD(5) 

A1,A2 

Ai,A2 

thermal  recovery  coefficient  for  hardening, 
(determined  by  CTD(6)  and  CTD(7)) 

AMI 

mi 

isotropic  hardening  rate  coefficient,  Cn(l) 

Z1 

Zl 

limiting  value  of  isotropic  hardening,  CTI(2) 

R1,R2 

th  ri 

thermal  recovery  exponent  for  hardening, 

cn(3) 

DO 

Do 

limiting  inelastic  strain  rate,  CTI(4) 

ANU 

V 

Poisson’s  ratio,  CTKS) 

G 

G 

shear  modulus 

AK 

K 

bulk  modulus 

DAM2 

laai 

differential  am2  with  solution  increment 

DAN 

differential  an  with  solution  increment 

DZO 

differential  Zq  with  solution  increment 

DZ2 

differential  Z2  with  solution  increment 

DZ3 

HDI 

differential  Z3  with  solution  increment 

DAI 

dA^  ^ 

a**' 

differential  Ai  with  solution  increment 

DA2 

ai  dt 

differential  A2  with  solution  increment 

DG 

dG  ar,^ 

differential  shear  modulus  with  solution 

increment 

DK 

aK  ai,^ 
ai  a'" 

differential  bulk  modulus  with  solution 

increment 
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Table  5 


State  Variables  of  Common  Blocks  BPSTAT2  and  BPSTAT3 


Variable  Name 

Symbol 

Variable  Description 

EIN 

wmm 

components  of  inelastic  strain 

ZI 

ZI 

current  isotropic  drag  stress 

SIGEFF 

V3J2 

effective  stress 

ZD 

ZD 

current  magnitude  of  directional  drag  stress 

BETA 

Pij 

components  of  directional  drag  stress 

EPSO 

— 

mechanical  strain  at  the  beginning  of 

ADINA's  time  TAU 

4.2  Internal  Variable  Names 

This  section  summarizes  the  internal  variables  of  the  user-defined 
subroutines  described  in  the  previous  section  (4.1).  Table  6  presents  ADINA6.0 
internal  variables  which  enter  into  CUSER2  and  CUSER3  from  subroutines 
USER2  and  USER3,  respectively.  Table  7  reviews  the  variables  found  within 
CUSER2  and  CUSER3,  while  Table  8  lists  the  variables  found  in  DBODNER2  and 
DBODNER3. 

Some  variables  change  values  at  certain  locations  within  the  numerical 
scheme.  For  example,  the  variable  STRESS  may  have  a  different  value  when 
exiting  subroutines  CUSER2  than  entering.  The  values  of  such  variables  depend 
on  solution  control  variables  --  TIME,  TAU,  INTER,  and  KTR.  The  variable  TIME 
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refers  to  the  time  at  the  beginning  of  a  major  time  step,  NSTEP.  TAU  is  the  value 
of  time  at  the  beginning  of  a  sub-incremental  time  step.  INTER  and  KTR  are 
control  variables  for  ADINA's  sub-incremental  solution  scheme.  ADINAIN3.0 
manuals  [17]  provide  details  of  this  sub-incremental  solution  scheme.  Variables, 
DT  and  DTAU,  refer  to  major  and  sub-incremental  time  steps  sizes,  respectively. 


Table  6 

ADIN  A  Variables  Supplied  to  CUSER2  and  CUSER3 


Variable 

Name 

Variable 

Description 

KEY  =  1 

initialize  the  working  arrays  during  input  phase 

KEY  =  2 

calculate  element  stresses 

KEY  =  3 

calculate  the  stress/strain  matrix 

KEY  =  4 

print  stresses  and  other  desired  variables  during  stress  print-out 

NG 

element  group  number 

NEL 

element  number 

IPT 

integration  point  number 

IT2D 

2d  element  type  identifier 

IT2D  EQ.O 

IT2D  EQ.l 

plane  strain  (in  y-z  plane) 

IT2DEQ.2 

plane  stress  (in  y-z  plane) 

IT2D  EQ.3 

plane  stress  (in  3/d  space) 

STRESS 

stress  components,  received  at  time  TAU  and  updated  to  time 
TAU+DTAU 

EPS 

total  strain  components  at  time  TIME. 

STRAIN 

total  strain  components  at  time  TIME+DT. 

DEPS 

subdivided  incremental  mechanical  strain  components  (total 
incremental  strain  minus  incremental  thermal  strains) 

DEPST 

components  of  sub  incremental  thermal  strain 

THSTRl 

total  thermal  strain  at  time  TAU 

THSTR2 

total  thermal  strain  at  time  TAU+DTAU 

INTER 

number  of  subdivisions  for  the  strain  increments 
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Table  6  (continued) 

ADINA  Variables  Supplied  to  CUSER2  and  CUSER3 


Variable 

Variable 

Name 

Description 

KTR 
KTR  =  1 
KTR  = 
INTER 

current  subdivision  number 
for  first  subdivision 
for  last  subdivision 

SCP 

solution  control  parameters 

SCP(1) 

relaxation  factor 

SCP(2) 

Bodner-Partom  convergence  tolerance 

SCP(3) 

element  number  for  special  history  output 

SCP(4) 

Gauss  point  for  special  history  output 

ARRAY 

working  array  (real)  for  user-defined  state  variables  at  time  TAU 
and  updated  by  user-supplied  coding  to  correspond  to  time 
TAU+DTAU 

lARRAY 

working  array  (integer)  for  user-defined  state  variables  at  time 
TAU  and  updated  by  user-supplied  coding  to  correspond  to  time 
TAU+DTAU 

D 

stress /strain  matrix ,  to  be  calculated  by  user-supplied  coding 

ALFA 

coefficient  of  thermal  expansion  at  time  TAU 

CTD 

temperature-dependent  material  constants  at  time  TAU 

CTD(l) 

elastic  modulus,  E 

CTD<2) 

kinetic  parameter,  n 

CTD(3) 

fully  recovered  (minimum)  value  of  isotropic  hardening,  Zo,Z2 

CTD(4) 

limiting  (maximum)  value  of  directional  hardening,  Z3 

CTD(5) 

hardening  rate  coefficient  of  directional  hardening,  m2 

ALFAA 

coefficient  of  thermal  expansion  at  time  TAU+DTAU 

CTDD 

temperature-dependent  material  constants  (same  as  CTD)  at  time 
TAU+DTAU 

cn 

temperature-independent  material  constants 

CTI(l) 

isotropic  hardening  rate  coefficient,  mi 

cn(2) 

limiting  value  of  isotropic  hardening,  Zi 

CTI(3) 

recovery  coefficient  for  drag  hardening,  ri,r2 
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Table  6  (continued) 

ADINA  Variables  Supplied  to  CUSER2  and  CUSER3 


Variable 

Variable 

Name 

Description 

CTI(4) 

limiting  inelastic  strain  rate.  Do 

cn(5) 

Poisson's  ratio,  v 

0X1(6) 

coefficient  for  recover  coefficient  (functional  form  for  Ai,  A2) 

CTI(7) 

exponential  for  recover  coefficient  (functional  form  for  A],  A2) 

TMPl 

TMP2 

temperature  at  integration  point  ITT  at  time  TAU+DTAU 

TIME 

time  at  current  step 

DT 

time  step  increment 

INTEG 

integration  parameter  for  stress  integration 

ISUBM 

flag  for  continuation  of  subdivision  in  time  step 

INDNL 

flag  for  element  formulation 
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Table? 

Internal  Variables  in  CUSER2  and  CUSER3 


Variable 

Variable 

Name 

Description 

RELAX 

relaxation  factor,  SCP(l) 

TOLER 

1  Bodner-Partom  convergence  tolerance,  SCP(2) 

PUNT 

Gauss  point  for  special  history  output,  SCP(3) 

HELE 

element  number  for  special  history  output,  SCP(4) 

E 

elastic  modulus,  CTD(l) 

AN 

kinetic  parameter,  CTD(2) 

Z0,Z2 

fully  recovered  (minimum)  value  of  isotropic  hardening,  CTD(3) 

Z3 

limiting  (maximum)  value  of  directional  hardening,  CTD(4) 

AM2 

hardening  rate  coefficient  of  directional  hardening,  CTD(5) 

AMI 

isotropic  hardening  rate  coefficient,  CTKl) 

Z1 

limiting  value  of  isotropic  hardening,  CTI(2) 

R1,R2 

recovery  exponent  for  hardening,  Cn(3) 

DO 

limiting  inelastic  strain  rate,  Cn(4) 

ANU 

Poisson's  ratio,  CTKS) 

DTAU 

time  inaement  step 

G 

shear  modulus 

AK 

bulk  modulus 

ALPHAl 

biaxial  model  parameter  (not  used) 

DAM2 

differential  am2  with  solution  increment 

DAN 

differential  an  with  solution  increment 

DZO 

differential  Zq  with  solution  increment 

DZ2 

differential  Z2  with  solution  increment 

DZ3 

differential  Z3  with  solution  increment 

DAI 

differential  Ai  with  solution  increment 

DA2 

differential  A2  with  solution  increment 

DG 

differential  shear  modulus  with  solution  increment 

DK 

differential  bulk  modulus  with  solution  increment 
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Table  8 

Internal  Variables  in  DBODNER2  and  DBODNER3 


Variable 

Variable 

Name 

Description 

DEVEPS 

deviatoric  strains 

EPS 

mechanical  strains  (total  minus  thermal) 

AVGEPS 

average  (mean)  normal  strain 

DDEPS 

sub-incremental  deviatoric  strains 

DEIEFF 

effective  inelastic  strain  increment 

DEPS 

incremental  mechanical  strain  (total  minus  thermal) 

DEPSI 

incremental  inelastic  deviatoric  strain 

DEVSIG 

deviatoric  stress 

ICOUNT 

iteration  loop  counter 

ISUB 

number  of  sub-increments  within  DBODNER2 

RELAX 

relaxation  factor  for  new  stress  estimate 

TOLER 

tolerance  for  convergence 

SIGNEW 

effective  stress  of  previously  converged  stress  state 

SIGOLD 

old  estimate  effective  stress 

SIGEST 

new  effective  stress  estimate 

SIGHYD 

hydrostatic  (mean)  stress 

STRESO 

stress  of  previously  converged  stress  state 

STRESS 

current  stress 

TIME 

current  time  value 

EINEST 

estimated  inelastic  strains 

EIN0(4) 

inelastic  strains  of  previously  converged  stress  state 

BETADOT 

directional  drag  stress  rate  vector 

BETAEST 

estimated  directional  drag  stress  vector 

BETA0(4) 

directional  drag  stress  of  previously  converged  stress  state 

u,v 

directional  unit  vectors 

ZIEST 

estimated  isotropic  drag  stress 
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4.3  Numerical  Algorithms 

The  algorithxns  that  solve  the  Bodner-Partom  constitutive  equations 
coiwist  of  a  mixture  of  iteration  loops  and  sub-incremental  schemes.  This  hybrid 
combination  of  iterations  and  sub-incrementation  works  well  for  the  inherently 
"stiff  nature  of  the  Bodner-Partom  differentials  equations  and  P21S  material 
parameters.  The  algorithms  discussed  in  this  section  are  found  in  subroutines 
DBODNER2and  DBODNER3. 

The  algorithm,  as  shown  in  Figure  1,  consists  of  two  iteration  loops  and  a 
sub-inaemental  solution  scheme.  Prior  to  sub-incremental  integration,  state 
variable  values  assume  their  pre-incremental  values;  the  number  of  sub¬ 
increments  (ISUB)  initializes  to  one;  and  constant  rate  variables  get  scale  by 
TF ACTOR.  The  integration  of  the  sub-incremental  loop  begins  at  Step  III.  The 
primary  iteration  loop  starts  at  step  III.D  and  converges  on  stress,  inelastic  strain, 
isotopic  hardening  and  the  directional  hardening  parameters.  Preliminary 
investigations  on  convergence  found  that  the  rate  equations  for  Pij,  Eq.  (6) 
strongly  depend  on  the  current  value  of  Pij,  thus  making  convergence  on  pij 
difficult.  The  value  for  often  diverges  during  the  primary  iteration  loop.  The 
second  iteration  loop,  starting  to  Step  III.D.3,  saves  computational  time  by 
determining  divergence  more  quickly  than  the  primary  iteration  loop,  since 
the  second  loop  avoids  several  numerical  operations  found  in  the  primary 
iteration  loop.  Non-convergent  solutions,  as  defined  by  a  maximum  limit  set  on 
iteration  steps,  return  to  the  beginning  of  the  sub-incremental  integration  (Step  I) 
with  an  increase  in  the  number  of  sub-increments,  ISUB.  When  the  maximum 
number  of  sub-incremental  step  equals  128,  the  solution  solver  ceases  operation 
and  prints  a  diagnostic  debugging  output. 
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I.  Store  variable  at  the  beginning  of  solution  increment 

A.  state  variables 

B.  temperature-dependent  material  parameters 

n.  Initialize  variables  for  sub-incremental  solution  step  cutting  of  KUB 

A.  initialize  estimated  new  values  for  state  variables 

B.  restore  material  parameters  from  Step  I 

C.  determine  new  material  parameters  rates  by  factors  of  ISUB 

D.  restore  state  variables  to  values  from  Step  I 

in  Begin  sub-incremental  solution  integration 

A.  update  all  temperature-dependent  parameters  to  end  of  sub- 
incremental  step 

B.  update  deviatoric  and  mean  mechanical  strains 

C.  step  iteration  counter  ICOUNT 

D.  begin  primary  iteration  loop 

1.  estimate  inelastic  strain  inaement  (engineering  inelastic 
shear 

strains  are  computed) 

2.  compute  new  stress  state  and  inelastic  work  rate 

3.  begin  secondary  iteration  loop  on  directional  hardening 
variable  BETA 

a.  compute  beta  rate 

b.  compute  new  estimate  for  beta 

c.  check  for  convergence  of  beta 

(1)  if  converged,  then  continue 

(2)  if  not,  then  increase  cutting  factor  ISUB  by 
factor  of  2.0  and  precede  to  Step  II.  A 

d.  stop  solution  for  excessive  number  of  time  cuts 

4.  compute  new  estimate  for  isotropic  hardening 

5.  check  for  convergence  of  effective  stress,  incremental 
inelastic 

strains,  isotropic  and  directional  hardening 

a.  if  converged,  then  continue  with  conclude  current  sub¬ 
increment  via.  STEP  III.E 

b.  if  not,  make  new  estimate  for  effective  stress  and 
increase  iteration  count  ICOUNT  BY  1. 

c.  check  for  excessive  iteration  count 

(1)  if  excessive,  then  increase  cutting  factor  ISUB 
by  factor  of  2.0  and  precede  to  Step  II. A 

(2)  if  not,  to  Step  ni.D 

E.  Update  converged  solution  with  estimates,  then  continue  a 
STEP  III.A 

IV.  Complete  all  sub-incremental  cycles  and  then  return  to  CUSER2  (or 
CUSER3) 


Figure  1  Numerical  Algorithm  for  Solving  of  Bodner-Partom  Equations. 
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Section  5 


Example.  Verification  and  Errors 

This  section  provides  an  example  of  the  application  of  the  Bodner-Partom 
user-defined  subroutines  with  the  P21S  material  parameters.  The  example 
ADINAIN3.0  input  file  for  a  simple  test  case.  This  section  also  presents  several 
numerical  exercises  that  serve  as  a  basis  for  comparisons  of  results  from  an 
independent  solution  source  with  those  from  the  finite  element  method.  The 
difference  between  the  results  from  these  two  solution  sources  are  measures  of 
accuracy  for  these  algorithms  and  assoiated  intergration  schemes.  The  last  part 
of  this  section  presents  an  investigation  that  guilds  the  users  of  these  subroutines 
to  minimize  errors  effectively. 

5.1  Numerical  Example 

The  file,  vcasel.in,  (see  Appendix  C)  is  an  example  ADINAIN3.0  input  file 
that  contains  the  required  P21S  material  parameters.  ADINA6.0  preprocessor 
ADINAIN3.0  generates  a  data  file  that  executes  the  user-defined  subroutines  that 
were  implemented  into  ADINA6.0.  The  example  consists  of  an  axisymmetric 
element  unidirectionally  loaded  at  a  constant  strain  rate  of  833.3xl0'^/s  and 
constant  temperature  of  25®C,  which  has  the  same  conditions  of  verification  test 
case  1  (discussed  below).  Input  of  the  P21S  material  properties  starts  with  the 
string  "MATERIAL  1  USER ..."  (see  Appendix  C).  The  cards  that  refer  to  material 
and  solution  parameters  are  described  in  Table  9  [17].  In  these  example  analyses, 
the  thermal  expansion  coefficients  are  set  to  zero,  since  the  comparisons  made 
here  only  consider  mechanical  strain  (total  strain  minus  thermal  strain). 
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Table  9  [17] 

Variable  Names  in  ADINAIN3.  0  Input  File 


Variable 

Name* 

Variable 

Description 

ALFA 

coefficient  of  thermal  expansion 

CTD(l) 

elastic  modulus,  E 

CTD(2) 

kinetic  parameter,  n 

CTD(3) 

fully  recovered  (minimum)  value  of  isotropic  hardening,  Zo,Z2 

CTD{4) 

limiting  (maximum)  value  of  directional  hardening,Z3 

CTD(5) 

hardening  rate  coefficient  of  directional  hardening,  m2 

CTI(1) 

isotropic  hardening  rate  coefficient,  mi 

CTI(2) 

limiting  value  of  isotropic  drag  stress,  Zi 

CTI(3) 

thermal  recovery  exponent  for  hardening,  ri,  r2 

CTI(4) 

limiting  inelastic  strain  rate,  Dq 

CTI(5) 

Poisson’s  ratio, v 

CTI(6) 

coefficient  for  thermal  recovery’s  functional  form  of  Ai  and  A2 

CTI(7) 

exponential  for  thermal  recovery’s  functional  form  for  Ai  and  A2 

XINTER 

number  of  incremental  subdivisions  of  NSTEP 

SCP(l) 

relaxation  factor 

SCP(2) 

Bodner  Partom  convergence  tolerance,  TOLER 

SCP(3) 

element  number  for  special  history  output  ’fort.53’ 

SCP(4) 

Gauss  point  for  special  history  output  file  ’fort.53' 

*User-defined  material  parameters  as  described  in  ADINAIN3.0  manual  [17] 
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This  section  presents  five  comparisons  of  the  finite  element  results  with 
those  obtained  through  an  independent  solution  source.  These  test  cases  consist 
of  three  isothermal  and  two  nonisothermal  simulations.  The  accuracy  measure 
for  this  investigation  is  percent  deviation  of  the  finite  element  stress  solution 
from  Mathematica  [16],  which  has  an  associated  error  of  approximately  10"^ 
percent.  The  percent  stress  deviation,  %D,  is  written  as: 

%D  =  100%  *  (9) 

<5math 

Where  CpE  and  Cmath  are  the  stresses  obtained  from  the  finite  element  and 
Mathematica  solutions,  respectively.  The  two-  and  thice-dimensional  finite 
element  stress  solutions  are  identical,  thus  only  one  finite  element  solution 
deviation  is  shown. 

All  test  cases  simulate  monotonically  loaded  elements  with  constant 
mechanical  strain  rates  of  either  833.3xl0‘6/s  or  8.33el0-6/s.  Table  10 
summarizes  all  the  test  conditions  and  location  in  the  appendices  of  the 
Mathematica  solution  files  used  in  the  comparisons.  The  data  in  these  files 
appear  in  Appendices  D  through  H.  Figures  2, 3  and  4  illustrate  the  stress-strain 
response  and  the  finite  element  stress  deviation  from  Mathematica  for  isothermal 
Test  Cases  1,  2,  and  3,  respectively.  Overall,  the  isothermal  load  cases  show  that 
the  finite  element  and  Mathematica  solutions  are  well  within  0.5  percent  of 
predicting  the  same  response. 


Table  10 


Verification  Test  Cases 


Test 

Case 

Number 

Temperature 

(X) 

Thermal 

Condition 

Strain  Rate 
(10-6/s) 

Mathermtica 
Solution  File 

Appendix 

1 

25 

isothermal 

833.3 

vcasel.math 

D 

2 

650 

isothermal 

833.3 

vcase2.math 

E 

3 

650 

isothermal 

8.333 

vcase3.math 

F 

4 

25/482/25 

non- 

isothermal 

833.3 

vcase4.math 

G 

5 

650/760/650 

non- 

isothermal 

833.3 

vcaseS.math 

H 

s. 

2 

w 

(0 

s 

03 


Mechanical  Strain 


Stress-Strain  Response  and  Finite  Element  Solution  Deviation  from 
Mathermtica  for  Test  Case  1 


Percent  Solution  Deviation 


Mechanical  Strain 


Figure  3  Stress-Strain  Response  and  Finite  Element  Solution  Deviation  from 
Mathermtica  for  Test  Case  2 
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Strain 

Figure  4  Stress-Strain  Response  and  Finite  Element  Solution  Deviation  from 
Mathermtica  for  Test  Case  3 
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The  nonisothermal  test  cases  are  constant  strain-rate  tensile  simulations 
with  changes  in  temperatures  at  prescribed  strain  levels.  These  non-isothermal 
loading  profiles  provide  an  excellent  basis  to  evaluate  the  various  forms  of  state- 
variable  temperature  differentials  (e.g.,  and  ^r)-  The  temperature 

profile  and  stress-strain  response  of  Test  Case  4  are  illustrated  in  Figure  5  with 
the  associated  finite  element  solution  deviation  from  Mathematica  presented  in 
Figure  6.  The  temperature  profile  and  stress  response  for  Test  Case  5  are  shown 
in  Figure  7  with  solution  deviation  given  in  Figure  8.  Overall,  the  finite  element 
and  mathematica  stresses  are  well  within  0.5  percent  of  predicting  the  same 
response  for  the  nonisothermal  segments  of  the  loading.  During  the 
nonisothermal  periods  of  loading  the  percent  deviation  increases  slightly,  which 
is  more  evident  in  Test  Case  5  than  in  of  Test  Case  4.  For  practical  use,  more 
solution  increments  and/or  smaller  convergence  tolerances  may  be  required  to 
achieve  the  same  level  accuracy  for  non  isothermal  loading  segments  as  those 
obtained  during  isothermal  periods,  especially  at  elevated  temperatures. 
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Mechanical  Strain 


Temperature  Profile  and  Stress  Response  for  Non-Isotherm 
Case  4 


Mechanical  Strain 


Stress-Strain  Response  and  Finite  Element  Solution  Deviation  from 
Mathematica  for  Test  Case  4 
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Temperature  Profile  and  Stress  Response  for  Non-Isothermal  Test 
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Figure  8  Stress-Strain  Response  and  Finite  Element  Solution  Deviation  from 
Mathermtica  for  Test  Case  5 
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Preliminary  investigations  with  Test  Case  5  revealed  an  unacceptable 
deviation  between  finite  element  and  Mathematica  solutions  (2.9%).  The  poor 
correlation  occurred  when  Test  Case  5  consisted  of  48  solution  increments,  rather 
than  the  96  increment  case  shown  in  Figure  8.  The  nature  of  the  excessive 
solution  deviation  is  determined  as  each  solution  parameter  is  systematically 
changed.  The  investigation  considers  two  types  of  solution  parameters  -- 
iteration  convergence  tolerances  and  solution  increment  number. 

Iteration  convergence  tolerances  occur  at  two  different  levels  within  the 
numerical  algorithms  --  global  and  local.  Global  tolerance  DNORM  is  a  measure 
of  convergence  on  displacements  and  forces  (energy)  within  the  entire  system  of 
nodes  and  elements  for  each  major  time  step,  NSTEP.  On  the  local  level,  within 
the  Bodner-Partom  iteration  schemes  (see  Section  4.3)  the  tolerance  parameter 
TOLER  measures  the  accuracy  of  the  converged  stress  state  for  strain  increments 
send  to  each  elemental  Gauss  point. 


The  global  parameter  DNORM  is  not  a  direct  measure  of  solution 
convergence  tolerance  within  ADENA.  Convergence  occurs  when  the  change  in 
displacement,  |du‘|,  for  each  global  iteration  estimate  becomes  relatively  small 
compared  to  DNORM,  which  is  expressed  as: 


d  u‘ 


DNORM 


<DTOL. 


Eq.  (10) 


With  the  displacement  tolerance  DTOL  set  constant,  a  change  in  DNORM 
directly  alters  the  solution  tolerance.  Since,  DNORM  parameter  has  no  standard 
default  value  within  ADINA,  this  investigation  considers  variations  in  DNORM 
rather  than  DTOL.  The  previously  completed  analyses  of  Section  5.2  used  a 
DNORM  base-line  value  of  l.OxlO"^.  Values  of  DNORM  for  this  investigation  are 
1.0xl0"3  and  1.0x10^.  The  change  of  DNORM  to  1.0x10  ^  or  1.0x10"^  produces  a 
only  a  small  change  in  stress  solution  (less  than  0.1%). 

The  changes  in  the  localized  tolerance  parameter  TOLER  within  the 
Bodner-Partom  iteration  loops,  results  in  minor  shift  in  the  solution  (0.25%),  as 
illustrated  in  Figure  9.  Overall,  the  tighter  tolerances  of  DNORM  and  TOLER 
still  produced  unacceptable  levels  of  solution  deviation  from  Mathermtica.  Note 
that  in  change  in  TOLER  from  0.001  to  0.00001  did  increase  the  finite  element 
solution  running  time  from  122  to  431  CPU  seconds*,  respectively,  as  shown  in 
Table  11.  Based  on  the  marked  debit  in  computational  efficiency  with  little 
increase  in  solution  accuracy,  increasing  TOLER  for  a  more  accurate  solution  is 
not  practical.  However,  the  TOLER  at  and  above  0.001  can  produce  non- 
convergent  solutions,  especially  if  the  finite  element  configuration  contains 
additional  non-linearities  (e.g.  gap  elements).  A  smaller  values  for  TOLER  to 
may  be  required  to  achieve  convergence  for  these  highlj'  non  linear  soh.’ticns. 


*  CPU  time  on  a  Sun  Sparc  Station  II  with  4/75  processor 
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Figvire  9  Effects  of  Localized  Tolerance  Parameter,  TOLER,  on  Finite 
Element  Solution  Deviation  from  Mathematica. 

ADINA's  integration  scheme  has  two  different  levels  of  incrementing  the 
solution  "  major  and  minor.  The  global  convergence  of  displacements  occurs  at 
the  major  solution  increment,  NSTEP.  The  minor  solution  increment  sub-divides 
the  major  increment  by  XINTER  for  quicker  integration  of  the  nonlinear  material 
models.  Increasing  either  NSTEP  or  XINTER  decreases  Ihe  solution  deviation,  as 
shown  in  Figures  10  and  11,  respectively.  The  decrease  of  solution  deviation 
with  increases  in  solution  time  steps  is  more  significant  than  changes  caused  by 
increasing  the  tolerance  parameters. 
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Figure  10  Effect  of  Major  Time  Step  NSTEP  on  Finite  Element  Solutic 
Deviation  from  Mathematica. 
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Figure  11  Effect  of  Sub-incremental  Parameter  XINTER  on  Finite  Element 
Solution  Deviation  from  Mathematica. 
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Percent  Solution  Deviation  3  Percent  Solution  Deviation 


To  summarize,  this  error  investigation  only  considers  a  simple 
unidirectional  test  case  during  nonisothermal  loading;  however,  knowledge  of 
the  results  will  save  time  and  effort  when  conducting  similar  error  investigation 
for  more  complex  cases.  First,  note  that  the  errors  associated  with  isothermal 
conditions  are  quite  small  compared  to  those  of  nonisothermal  loading  periods, 
even  with  a  few  time  steps.  For  the  nonisothermal  Test  Case  5,  Table  11  presents 
a  summary  of  finite  element  solution  parameters  and  their  effects  on  computer 
CPU*  time  and  solution  deviation.  Increases  in  solution  tolerances  DNORM  and 
TOLER  do  not  significantly  change  the  resulting  stress  error;  however,  CPU  time 
did  increase.  Increases  in  number  of  solution  increments,  NSTEP  or  XINTER, 
significantly  decrease  the  solution  error,  while  CPU  time  remains  relatively 
constant.  Overall,  the  most  accurate  solution  per  CPU  time  occurs  with  increased 
number  of  major  time  steps,  NSTEP. 


*CPU  time  on  a  Gun  Sparc  Station  H  with  4/75  processor 
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Table  11 


Effects  of  Finite  Element  Solution  Parameters 
on  CPU*  Time  and  Solution  Deviation 


NSTEP 

DNORM 

XINTER 

TOLER 

CPU 
Time  (s) 

Maximum 

Percent 

Deviation 

48 

1.0x10-5 

10 

0.0001 

206 

2.96 

48 

1.0x10-4 

10 

0.0001 

173 

3.02 

48 

1.0x10-5 

10 

0.0001 

149 

3.02 

48 

1.0x10-5 

5 

0.0001 

180 

6.51 

48 

1.0x10-5 

20 

0.0001 

204 

1.25 

48 

1.0x10-5 

10 

0.001 

122 

3.16 

48 

1.0x10-5 

10 

0.00001 

432 

3.23 

96 

1.0x10-5 

10 

0.0001 

200 

1.25 

192 

1.0x10-5 

10 

0.0001 

179 

0.362 

384 

1.0x10-5 

10 

0.0001 

231 

0.0478 

Mathermtica 

— 

— 

— 

2638 

— 

CPU  time  on  a  Sun  Sparc  Station  I[  with  4/75  processor 


35 


Section  6 


References 


(1)  Bodner,  S.R.,  and  Partom,Y.,  "Constitutive  Equations  for  Elastic- 
Viscoplastic  Strain  Hardening  Materials,"  ASME  Journal  of  Applied 
Mechanics,  Vol.  42,1975,  p.385. 

(2)  Bathe,  KJ.,  "Finite  Elements  in  CAD  and  ADINA,"  Nuclear  Engineering 
and  Design,  Vol.  98, 1986,  pp.57-67. 

(3)  Chan,  K.S.,  Bodner,  S.R.,  and  Lindholm,  U.S.,  "Phenomenological 
Modeling  of  Hardening  and  Thermal  Recovery  in  Metals,"  ASME  Journal 
of  Engineering  Materials  and  Technology,  Vol.  110,  January  1988,  pp.1-8. 

(4)  Chan,  K.S.,  Lindholm,  U.S.,  Bodner,  S.R.,  and  Nagy,  A.,"High  Temperature 
Inelastic  Deformation  of  the  B1900+Hf  Alloy  Under  Multiaxial  Loading; 
Theory  and  Experiment,"  ASME  Journal  of  Engineering  Materials  and 
Technology,  Vol.  112,  January  1990,  pp.  7-14. 

(5)  Chan,  K.S.,  Lindholm,  U.S.,  Bodner,  S.R.,  "Constitutive  Modeling  for 
Isotropic  Materials  (HOST),  NASA  CR-182132  (SwRI-7576),  NASA  Lewis 
June  1988. 

(6)  Chan,  K.S.  and  Lindholm,  U.S.  "Inelastic  Deformation  Under 
Nonisothermal  Loading,"  ASME  Journal  of  Engineering  Materials  and 
Technology,  Vol.  112,  January  19k),  pp.15-25. 

(7)  Sherwood,  J.A.  and  Boyle,  M.J.,  "Investigation  of  the  Thermomechanical 
Response  of  a  Titanium  Aluminide /Silicon  Carbide  Composite  Using  a 
Unified  State  Variable  Model  in  ADINA,"  Computers  and  Structures,  Vol. 
40, 1991,  pp.257-269. 

(8)  Ramaswamy,  V.G.,  Stouffer,  D.C.,  and  Laflen,  J.H.,  "A  Unified 
Constitutive  Model  for  the  Inelastic  Uniaxial  Response  of  Rene'80  at 
Temperatures  Between  538C  and  982C,"  Transactions  of  the  ASME,  Vol. 
112,  July  1990,  pp.  280-286. 

(9)  Stouffer,  D.C.,  Ramaswamy,  V.G.,  Laflen,  J.H.,  Van  Stone,R.H.,  and 
Williams,  R.,  "A  Constitutive  Model  for  the  Inelastic  Multiaxial  Response 
of  Rene’80  at  871C  and  982C,"  ASME  Journal  of  Engineering  Materials  and 
Technology,  Vol.  112,  April  1990,  pp.241-246. 


36 


References  (Continued) 


(10)  McDowell,  D.L.,  "A  Bounding  Surface  Theory  for  Cyclic 
Thermoplastidty,"  ASME  Journal  of  Engineering  Materials  and 
Technology,  Vol.  114,  July  1992,  pp.  297-303. 

(11)  Walker,  K.  P.,  "Research  and  Developement  Program  for  Nonlinear 
Structuaral  Modeling  with  Advanced  Time-Temperature  Dependent 
Comstitutive  Relationships,"  NASA  Report  CR-165533, 1981. 

(12)  Freed,  A.  D.,  "Thermoviscoplastic  Model  With  Application  to  Copper," 
NASA  Report  TP-2845, 1988. 

(13)  Chaboche,  J.  L.,  "Constitutive  Equations  for  Cyclic  Plasticity  and  Cyclic 
Viscoplastidty,"  International  Journal  of  Plastidty,  Vol.  5, 1989,  pp.  247- 
302. 

(14)  Wang,  J.-D.  and  Ohno,  N.,  'Two  Equivalent  Forms  of  Nonlinear  Kinematic 
Hardening:  Application  to  Nonisothermal  Plasticity,"  International 
Journal  of  Plastidty,  Vol.  7, 1991,  pp.  637-650. 

(15)  Ashbaugh,  N.E.,  and  Khobaib,  M.,  Unpublished  Data,  Universty  of 
Dayton  Research  Institute,  Dayton,  Ohio. 

(16)  Wolfram,  S.,  Mathematica-A  System  for  Doing  Mathematics  bv  Computer. 
Sec.  Ed.,  1991,  Addison-Westley  Publishing  Company,  Inc.,  Redwood 
City,  C  A. 

(17)  ADINA-IN  for  ADINA  Users  Manual.  Report  ARD  90-4,  ADINA  R&D. 
Inc.,  September  1990,  71  Elton  Ave,  Watertown,  MA.,  02172. 


37 


Appjendix  A 


User-Defined  Subroutines  for  Two-Dimensional 
Axisymmetric  Elements  (file  cuser2_dbeta.O 


A 

B 

C 

C*I 

C*I 

C*I 

C*I 

ei 

ei 

c*i 

c*i 

c*i 

c*i 

ei 

oi 

c*i 

ei 

c*i 

ei 

ei 

C*! 

c*i 

c*i 

c*i 

ei 

ei 

c*i 

ei 

ei 

ei 

ei 

ei 

c*i 

CT 

ei 

c*i 

c*i 

ei 

ei 

ei 

c*i 

c*i 

c*i 

ei 


SUBROUTINE  CUSER2  (NG,NEmPTJT2D3TRESS,EPS,STRAIN,DEPS,DEPST, 
THSTR1,THSTR2,KTRJNTER,SCP,ARRAY,IARRAY,D, 
ALFA,CTD,ALFAA,CTDD,CTI.TMP1,TMP2,TIME,DTAU, 
PHIST,PRST,ANGLE,DPSP,INTEG,ISUBM,1NDNL,KEY) 


THIS  SUBROUTINE  IS  TO  BE  SUPPLIED  BY  THE  USER  TO  CALCULATE 
THE  STRESSES  AND  CONSTITUTIVE  MATRIX  OF  A  SPECIAL  MATERIAL. 

THIS  SUBROUTINE  IS  CALLED  IN  USER2  FOR  EACH  INTEGRATION  POINT 
FOR  EACH  2-D  SOLID  ELEMENT  TO  PERFORM  THE  FOLLOWING  OPERATIONS  : 

KEY.EQ.1  INITIALIZE  THE  WORKING  ARRAYS  DURING  INPUT  PHASE 

KEY.EQ.2  CALCULATE  ELEMENT  STRESSES 

KEY.EQ.3  CALCULATE  THE  STRESS/STRAIN  MATRIX 

KEY.EQ.4  PRINT  CALCULATED  STRESSES  AND  OTHER  DESIRED 
VARIABLES  DURING  STRESS  PRINTOUT 


THE  FOLLOWING  VARIABLES  ARE  USED  TO  PERFORM  THE  ABOVE  OPERATIONS 

NG  ELEMENT  GROUP  NUMBER 

NEL  ELEMENT  NUMBER 

IPT  INTEGRATION  POINT  NUMBER 

rr2D  2D  ELEMENT  TYPE  IDENTIFIER  ( NPAR(5) ) 

EQ.O  AXISYMMETRIC  (IN  Y-Z  PLANE) 

EQ.l  PLANE  STRAIN  (IN  Y-Z  PLANE) 

EQ.2  PLANE  STRESS  (IN  Y-Z  PLANE) 

EQ.3  PLANE  STRESS  (IN  3/  D  SPACE) 

STRESS(4)  STRESS  COMPONENTS,  RECEIVED  AT  TIME  TAU 
AND  UPDATED  BY  USER-SUPPLIED  CODING  TO 
CORRESPOND  TO  TIME  TAU+DTAU 

EPS(4)  TOTAL  STRAIN  COMPONENTS  AT  TIMET. 

IN  CASE  OF  LARGE  STRAIN  FORMULATION 
(U.L.H.),  EPS(4)  ARE  ELASTIC  STRAINS  AT 
TIMET. 
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0 


ei 

ei 

ei 

ai 

ei 

ei 

ei 

ei 

ei 

ei 

ei 

ei 

c*i 

c*i 

CM 

c*i 

ei 

ei 

c*i 

ei 

ei 

ei 

ei 

ei 

c*i 

ei 

ei 

c*i 

ei 

ei 

ei 

CM 

C*! 

c*i 

c*i 

cr^i 

ei 

ei 

c*i 

CM 

CM 

CM 

CM 

CM 

C*I 

ei 

ei 

ei 

CM 

CM 

CM 

CM 


STRAIN(4)  TOTAL  STRAIN  COMPONENTS  AT  TIME  T+DT. 
INCASEOFU.L.H.: 

A)  FOR  KEY.EQ.2  -  ELASTIC  STRAINS 
AT  TIME  TAU+DTAU  (FOR  KTR.EQ.0  - 
TRIAL  ELASTIC  STRAINS) 

B)  FOR  KEY.GT.2  -  ELASTIC  STRAINS  AT 
TIME  T+DT 

DEPS{4)  SUBDIVIDED  INCREMENTAL  STRAIN  COMPONENTS 
DEPS(1)  =  ( STRAIN(I)  -  EPSd)  )/INTER 
-  DEPST(I) 

( PASSED  TO  SUBROUTINE  CUSER2  BY  THE 
PROGRAM  ADINA ) 

DPSP(4)  INCREMENT  OF  INELASTIC  STRAIN  ( PLASTIC 
AND/OR  CREEP  AND/OR  VISCOPLASTIC,  ETC.) 

IN  THE  SUBDIVISION.  CALCULATED  BY  USER- 
SUPPLIED  CODING  AND  PASSED  TO  THE  PRCX^RAM 
ADINA.  USED  FOR  U.L.H.  ONLY. 

DEPST(4)  COMPONENTS  OF  SUBINCREMENTAL  THERM  AL 

STRAIN 

( PASSED  TO  SUBROUTINE  CUSER2  BY  THE 
PROGRAM  ADINA ) 

THSTRl  (4)  TOTAL  THERMAL  STRAIN  AT  TIME  TAU 

FOR  KEY.EQ.4  -  THERMAL  STRAIN  AT  TIME  T 

THSTR2(4)  TOTAL  THERMAL  STRAIN  AT  TIME  TAU+DTAU 

FOR  KEY  EQ.4  -  THERMAL  STRAIN  AT  TIME  T+DT 

INTER  NUMBER  OF  SUBDIVISIONS  FOR  THE  STRAIN 

INCREMENTS  (INTER  =  INT(PROP(l  23)) 

KTR  CURRENT  SUBDIVISION  NUMBER 

EQ.0  CALCULATION  OF  TRIAL  ELASTIC 
STATE,  IN  CASE  INTEG*1 
EQ.l  FOR  HRST  SUBDIVISION 

EQ.INTER  FOR  LAST  SUBDIVISION 

SCP(4)  SOLUTION  CONTROL  PARAMETERS 

ARRAY(60)  WORKING  ARRAY  (REAL) ,  RECEIVED  AT 
TIME  TAU  AND  UPDATED  BY  USER-SUPPLIED 
CODING  TO  CORRESPOND  TO  TIME  TAU+DTAU 

IARRAY(2)  WORKING  ARRAY  (INTEGER) ,  RECEIVED  AT 

TIME  TAU  AND  UPDATED  BY  USER-SUPPLIED 
CODING  TO  CORRESPOND  TO  TIME  TAU+DTAU 

D(4,4)  STRESS/STRAIN  MATRIX ,  TO  BE  CALCULATED 


j 
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ei 

BY  USER-SUPPLIED  CODING 

ei 

C*I 

ALFA 

COEFFICIENT  OF  THERMAL  EXPANSION  AT 

C*I 

TIME  TAU 

C*I 

C*I 

CTD(5) 

TEMPERATURE-DEPENDENT  MATERIAL  CONSTANTS 

ei 

AT  TIME  TAU 

C*I 

C*I 

ALFAA 

COEFFICIENT  OF  THERMAL  EXPANSION  AT 

ei 

TIMETAU+DTAU 

C*I 

C*I 

(MDCKS) 

TEMPERATURE-DEPENDENT  MATERIAL  CONSTANTS 

C*I 

AT  TIMETAU+DTAU 

C*! 

C*I 

Cn(8) 

TEMPERATURE-INDEPENDENT  MATERIAL  CONSTANTS 

C*I 

ei 

TMPl 

TEMPERATURE  AT  INTEGRATION  POINT  IPT  AT 

ai 

TIME  TAU.  FOR  KEY.EQ.4  -  TEMPERATURE  AT 

C*I 

TIMET 

ei 

C*1 

TMP2 

TEMPERATURE  AT  INTEGRATION  POINT  IPT  AT 

C*I 

TIME  TAU+DTAU.  FOR  KEY.EQ.4  -  TEMPERATURE 

ei 

ATTIMET+DT 

C*I 

ei 

TIME 

TIME  AT  CURRENT  STEP ,  T+DT 

C*I 

ei 

DT 

TIME  STEP  INCREMENT ,  DT 

ei 

C*I 

PHIST(33) 

MATRIX  CONTAINING  DIRECTTION  COSINES  OF 

C*I 

PRINCIPAL  STRETCH  DIRECTIONS,  IN  CASE  OF 

C*I 

U.L.H.  FORMULATION 

C*I 

C*I 

PRST(3> 

PRINCIPAL  STRETCHES,  IN  CASE  OF  U.L.H. 

C*I 

ei 

ANGLE(2) 

ANGLES  OF  THE  FIRST  AND  SECOND  (IN  PLANE) 

ei 

PRINCIPAL  STRETCH  DIRECTIONS  WITH  RESPECT 

C*I 

TO  Y-AXIS,  U.L.H. 

CM 

C*I 

INTEG 

INTEGRATION  PARAMETER  FOR  STRESS 

C*1 

INTEGRATION 

(?»i 

EQ.0  -  FORWARD  INTEGRATION 

C*I 

EQ.I  -  BACKWARD  INTEGRATION 

C*I 

(RETURN  MAPPING) 

ei 

C*I 

ISUBM 

FLAG  FOR  CONTINUATION  OF  SUBDIVISION 

C*I 

IN  TIME  STEP,  APPLICABLE  FOR  INTEG.EQ.l 

C*I 

EQ.0  -  CONTINUATION 

CM 

EQ.-I  -  STOP  OF  SUBDIVISION 

C*1 

IN  THE  USER-SUPPLIED  CODING  THE  FLAG 

C*1 

(INITIALLY  EQ.0)  MUST  BE  SET  TO  -1 

ei 

WHEN  CRITERIA  FOR  STOPPING  SUBDIVISIONS 

ei 

ARE  REACHED 
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C*I 

C*I  INDNL  FLAG  FOR  ELEMENT  FORMULATION  (NPAR(3)) 

C*I  EQ.l  -  MATERIALLY  NONLINEAR  ONLY  (M.NO.) 

C*!  EQ.2- LARGE  DISPLACEMENTS  AND  SMALL 

C*I  STRAINS  (T.L.) 

ei  EQ.3- LARGE  DISPLACEMENTS  AND  LARGE 

C*I  STRAINS  (U.L.H.) 

C*I 

C*I 

C*I  NOTIE  THAT  THE  VARIABLES  PASSED  TO  THE  SUBROUTINE  WHEN  KEY=3,4 
OI  ARE  THESE  CALCULATED  LAST:  I.E.,  CALCULATED  IN  THE  LAST 
C*I  SUBDIVISION  WHEN  KEY=2.  HENCE  THESE  VARIABLES  ARE  NOT 
C*I  CALCULATED  WHEN  KEY=3,4. 

C*I 

C*I 

C*I 

IMPLICIT  DOUBLE  PRECISION  (  A-H,0-Z  ) 

C 

DIMENSION  STRESS(4)^TRAIN(4),DEPS(4).ARRAY(60),IARRAY(2),D(4,4), 

A  CTD(5),CTI(8),EPS(4).SCP(4) 

DIMENSION  CTDD(5),DEPST(4),THSTR1(4),THSTR2{4) 

DIMENSION  PHIST(3^LPRST(3),ANGLE(2).DPSP(4) 

C 

C 

REAL’S  DTAU 

C  TIME  INCREMENT  STEP,  DT 

REAL’S  E 

C  ELASTIC  MODULUS 

REAL’S  EIDEFF 

C  EFFECTIVE  INELASTIC  STRAIN  RATE 

REAL’S  G 

C  SHEAR  MODULOUS 

INTEGER  HINT 

C  INTEGRATION  POINT  WRITTEN  TO  HISTORY  FILE 

INTEGER  HELE 

C  ELEMENT  NUMBER  ASSOCIATED  WITH  HIPT 

INTEGER  I 

C  DO  LOOP  COUNTER 

INTEGER  J 

C  DO  LOOP  COUNTER 

C 

C  DIRECTIONAL  BODNER-PARTOM  MATERIAL  CONSTANTS: 

C 

C  TEMPERATURE  INDEPENDENT  CONSTANTS 
C 

REAL’S  AMI 
REAL’S  Z1 
REAL’S  R1 
REAL’S  R2 
REAL’S  DO 
REAL’S  ALPHAl 
C 
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C  TEMPERATURE  INDEPENDENT  CONSTANTS 
C 

REAL’S  AN 
REAL’S  AM2 
REAL’S  ZO 
REAL’S  Z2 
REAL’S  Z3 
REAL’S  A1 
REAL’S  A2 
C 

C  TEMPERATURE  DIFFERENTIALS 
C 

REAL’S  DAM2 
REAL’S  DAN 
REAL’S  DZO 
REAL’S  DZ2 
REAL’S  DZ3 
REAL’S  DAI 
REAL’S  DA2 
C 

COMMON  /MATCONST/  E,  G,  AK.  DG,  DK, 

&  AMI,  AM2,  ALPHAL  DAM2, 

&  Zl,  Z3,  Rl,  R2,  DO, 

&  AN,  ZO,  Zl,  Al,  A2, 

&  DAN,  DZO,  DZ2,  DZ3,  DAI,  DA2 

C 

C  COMMON  BLOCK  FOR  STATE  VARIBLES 
C 

COMMON  /BPSTAT2/  EIN(4),  Zl,  SIGEFF,  BETAEFF,  EPEFF,  ALPHA, 

&  ZD,BETA(4),EPS0(4) 

C 

C  RECALL  STATE  VARIBLES  AT  TIME  TAU 
C 

CALL  STGET2(ARRAY) 

C 

GO  TO  (1,23,4),  KEY 
ei 
O’! 

C*l  KEY  =  1 
CM 

<yi  INITIALIZE  COMPONENTS  OF  REAL  AND  INTEGER  WORKING  ARRAYS 
C’l  ( INITIALIZE  ARRAY(60)  AND  IARRAY(2) ) 

C*I 

I  CONTINUE 
C’l 

C’l  ”•  INSERT  USER-SUPPLIED  CODING 

C 

C 

C  INITIALIZE  STATE  VARIBLES 
C 

EIN(l)  =0.0 
EIN(2)  =0.0 
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EIN(3)  =0.0 
EIN(4)  =  0.0 
BETA(1)*0.0 
BETA(2>  =  0.0 
BETA(3)  =  0.0 
BETA(4)  =  0.0 
ZI  « CTD(3) 

SIGEFF  =  0.0 
EPEFF  =0.0 
BETAEFF  =  0.0 
ALPHA  =00 
ZD  =0.0 
C 

C  STORE  STATE  VARIBLES 
C 

CALL  STPUT2(ARRAY) 

C*I 

RETURN 

C*I 

OI 

C*I  KEY  =  2 

ei 

OI  INTEGRATION  OF  ELEMENT  STRESSES 
C*I  ( CALCULATE  STRESS(4) ) 

C*I 

2  CONTINUE 
C*I 

C*I  INSERT  USER-SUPPLIED  CODING 
C*I 

C  BODNER-P ARTOM  MATERIAL 
C 

C  FOR  2-D  SOLID  ELEMENT 
C 

C  GET  CURRENT  VALUES  OF  MATERIAL  CONSTANTS  FOR  GIVEN  TEMPERATURE 
C 

AMI  =cn(l> 

Zl  =  CTI(2) 

R1  =  cn(3) 

R2  =  CTI(3) 

DO  =  CTI(4) 

ANU  =CTI{5) 

ALPHAl  =  CTI(8) 

C 

At  =  Cn(6)  *  EXP(-Cn(7)/(TMP1  +273.)) 

A2  =A1 

DAI  =  CTI(6)  *  EXP{-CTI(7)/(TMP2+273.))-Al 
DA2  =  DAI 
C 

G  =  CTD0)/{2.*(1.+ ANU)) 

DG  =  CTDD(1)/(2.‘{1.+ ANU))-G 
C 

AK  =  CTD<l)/(l.-2.*  ANU) 
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DK  =  CTDD(l)/(l.-2  *  ANU)  -  AK 
C 

AN  «  CTD(2) 

DAN  =CTDD(2)-CTD(2) 

C 

22  =CTD(3) 

DZO  aCTDD(3)-CTD(3) 

DZ2  =  CTDD(3)-CTD(3) 

C 

Li  =  CTD(4) 

DZ3  »CTDD(4)-CrD(4) 

C 

AM2»  CTD(5) 

DAM2  =  CTDD(5)-CrD{5) 

C 

RELAX  =  SCP(1) 

TOLER  =  SCP(2) 

HINT  =  1FIX(SCP(3» 

HEUE  =  inX(SCP(4)) 

C 

C  COMPUTE  INCREMENTAL  STRAINS  WITH  INCREMENTAL  THERMAL  STRAINS 
C 

IF  (KTR  .EQ.  1)  THEN 
DO  20 1  =  1,4 

20  EPS0(I)  *  EPS(I)  -  THSTRKI) 

ENDIF 

C 

DTAUl  ®  DTAU  /  ROAT(  INTER ) 

C 

C  CALCULATE  INCREMENTAL  STRESSES  AND  UPDATE  STRESS  VECTOR 
C 

CALL  DBODNER2  ( DEPS,  STRESS,  DTAUl, 

+  NEL,IPT,  KTR,  TIME,  INTER,  RELAX,  TOLER) 

C 

DO  30  1  =  1,4 

30  EPS0(I)  =  EPSOd)  +  DEPS(I) 

C 

C - 

C  STORE  STATE  VARIBLES  AT  TIME  TAU  +  DTAU 
C 

CALLSTPUT2(ARRAY) 

C 

RETURN 

C»I 

C*I 

C*I 

C*I  KEY  =  3 
ei 

ei  FORM  CONSTITUTIVE  LAW 
cn  ( CALCULATE  D<4,4) ) 

CM 

3  CONTINUE 
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INSERT  USER-SUPPLIED  CODING 


C*1 
C*I 
C 

C- 
C 

C  UPDATE  THE  CONSTITUTIVE  MATRIX 
C 

c - 

c 

c - 

c 

C  THE  TEMPERATURE  DEPENDENT  ELASTIC  MATRIX  IS  DETERMINED 
C  WITHOUT  A  PLASTIC  CORRECTION 
C 

E  =CTDD(1) 

ANU  =  CTKS) 

C 

A1  =  E/(1.+ANU) 

Cl  =  Al*0.5 

A1  »  A1/(1.-2*ANU) 

B1  =  ArANU 
C 

DO  315 1=1,4 
D0315J=1,4 
315  D(I,J)  =  0 
C 

EX1,1)  =  A1 
IXU)  =  B1 
D(l,4)  =  B1 
D(2,l)  =  B1 
D(2,2)  =  A1 
D(2,4)  =  B1 
D(33)  =  Cl 
D(4,1)  =  B1 
D(4,2)  =  B1 
D(4,4)  =  A1 
C 

C  STORE  STATE  VARIBLES  AT  TIME  TAU  +  OTAU 
C 

CALL  STPUT2(ARRAY) 

C 

RE1URN 

C 

C*I 

C*I  KEY  =  4 
CM 

C*I  PRINTING  OF  ELEMENT  RESPONSE 
C*I  ( PRINT  STRESS(4),STRAIN(4) ) 

C*I 

4  CONTINUE 
CM 

CM  •**  INSERT  USER-SUPPLIED  CODING 
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ei 

ei 

c 

IF((NEL.EQ.HELE).AND.(IPT.EQ.HINT))THEN 
7006  FORMAT{2(lX,F8.1),]2(lX.E11.4)) 

ZMECHS  =  STRAIN(2)  -  THSTR2{2) 

ZTOT  =ZI+ZD 

WRrrE(53,7006)  time,tmp2,stress,einz:mechs, 

&  sigeff 

EKDIF 

C 

C  PRINT  HEADING  AND  ELEMENT  NUMBER 
C 

IF(IPT.EQ.l)  WRITE(6,9001)  NEL 
C 

C  PRINT  STRESSES 
C 

MODE  =  *  ELASTIC 

IF  (EIDEFF.GT.l  .E-9)  MODE='INELASTIC 
WRITE(6,9002)IPT3TRESS,SIGEFF,MODE3TRAIN,EIDEFF 
C 

C  FORMAT  STATEMENTS 
C 

9001  FORMAT  (//,4X,3HNEL,4X^HIPT,6X,9HSTRESS-YY,6X, 

+  9HSTRESS-ZZ,6X,9HSTRESS-YZ, 

+  6X,9HSTRESS-XX,4X,1 1 HEFF.  STRESS,/. 

+  20X,9HSTRAIN-yY,6X, 

+  9HSTRAIN-ZZ,6X,9HSTRAIN-YZ, 

+  6X,9HSTRAIN-XX,4X,6HEIDEFF,//,I7) 

9002  FORMAT  (7X,I7,5(2X,1PE13.6),3X,A9./,14X,5(2X,E13.6)) 
IF((NEL.EQ.1).AND.(NG.EQ.HNG).AND.{IPT.EQ.1))THEN 

REWINCKSS) 

:  NDIF 
C 

CALL  STPUT2(ARRAY) 

C 

RETURN 

END 

C 

SUBROUTINE  DBODNER2  (  DEPS,  STRESS,  DTAU, 

+  NEL,  IPT,  KTR,  TIME,  INTER,  RELAX,  TOLER) 

C 

IMPLICIT  DOUBLE  PRECISION  (  A-H,0-Z ) 

C - 

c 

REAL*8  DEVEPS{4) 

REAL'S  EPS(4) 

C 

REAL'S  AVGEPS 

C  AVERAGE  NORMAL  STRAIN  INCREMENT  AT  TIME  T 

REAL'S  DDEPS(4) 

C  SUBINCREMENTIAL  DEVITORIC  STRAINS 
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REAL*8  DEIEFF 

C  EFFECTIVE  INELASTIC  STRAIN  INCREMENT 

REAL’S  DEPS(4) 

C  INCREMENTAL  TOTAL  STRAIN  VECTOR 

REAL’S  DEPSI(4) 

INELASTIC  PORTION  OF  THE  DEVIATORIC  STRAIN 
VECTOR 

REAL’S  DEVSIG(4) 

DEVIATORIC  STRESS  VECTOR  OF  THE  STRESS  STATE 
AT  TIMET 
INTEGER  ICOUNT 

ITERATION  LOOP  COUNTER 
REAL’S  FACri 

INTEGER  FACTOR 
REAL’S  FACT2 

INTEGER  FACTOR 
INTEGER  INTER 

NUMBER  OF  SUBINCREMENTS 
INTEGER  IPT 

GAUSS  POINT 
INTEGER  ISUB 

NUMBER  OF  SUBINCREMENTS  WITHIN  C 
INTEGER  KTR 

CURRENT  SUBDIVISION  NUMBER 
EQ.  1  FOR  HRST  SUBDIVISION 
EQ.  2  FOR  LAST  SUBDIVISION 
INTEGER  NEL 

ELEMENT  NUMBER 
REAL’S  RELAX 

RELAXATION  FACTOR  FOR  NEW  STRESS  ESTIMATE 
REAL’S  TOLER 

TOLERANCE  FOR  CONVERGENCE 
REAL’S  SIGNEW 

NEW  EFFECTIVE  STRESS 
REAL’S  SIGOLD 

INITIAL  EFFECTIVE  STRESS 
REAL’S  SIGEST 

EFFECTIVE  STRESS  FOR  TIME  STEP  SIZE  OF  DT/JINT 
REAL’S  SIGHYD 

HYDROSTATIC  STRESS 
REAL’S  SIGMAX 

MAXIMUM  ABSOLUTE  VALUE  OF  EITHER  SIGEFF  OR  SIGEFZ 
FOR  INVESTIGATING  CONVERGENCE 
REAL’S  STRES0(4) 

C  STRESS  VECTOR  AT  TIME  TAU 

REAL’S  STRESS(4) 

C  STRESS  VECTOR  AT  TIME  TAU+DTAUT 

REAL’S  TIME 

C  CURRENT  TIME  VALUE 

REAL’S  EINEST(4) 

C  ESTIMATED  PLASTIC  STRAINS 

REAL’S  E1N0<4) 
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c 

PLASTIC  STRAINS  AT  TIME  TAU 

REAL’S 

BETADOT(4) 

c 

BETA  RATE  VECTOR 

REAL’S 

BETAEST(4) 

c 

ESTIMATED  BETA  VECTOR 

REAL’S 

BETA0(4) 

c 

BETA  VECTOR  AT  TIME  TAU 

REAL’S 

U(4),  V(4) 

c 

TEMPERATURE  VECTOR 

REAL’S 

ALEST 

c 

ESTIMATED  ALPHA 

REAL’S 

ZIEST 

c 

ESTIMATED  ISOTROPIC  DRAG  STRESS 

c 

c - 

c 

COMMON  /MATCONST/  E,  G,  AK,  DG,  DK, 

&  AMI,  AM2,  ALPHAl,  DAM2, 

&  Zl,  23,  Rl,  R2,  DO, 

&  AN,  ZO,  Z2,  Al,  A2, 

&  DAN,  DZO,  DZ2,  DZ3,  DAI,  DA2 

C 

COMMON  /BPSTAT2/  EIN(4),  ZI,  SIGEFF,  BETAEFF,  EPEFF,  ALPHA, 
&  ZD,  BETA(4),  EPS0{4) 

C 

SQRT3  «  SQRT(3.) 

C 

C  STORE  OLD  STATE  V  ARIBLES 
C 

DO  10 1=1,4 
STRESO(I)  =  STRESS(l) 

EIN0(I)  =EIN(I) 

BETAO(I)  =BETA(1) 

10  CONTINUE 
C 

ZIO  =ZI 
ALPHAO  =  ALPHA 
ZDO  =ZD 
C 

C  STORE  OLD  MATERIAL  CONSTANTS 
C 

GO  =  G 
AKO  =  AK 
ZOO  =  ZO 
Z20  =  Z2 
Z30  =  Z3 
AlO  =  Al 
A20  =  A2 
AM20=  AM2 
C 

C  INITIALIZE  OTHER  VARIBLES 
C 
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ISUB  =  1 
IDBUG  =0 

INITIALIZE  VARIABLES  FOR  SUB-TIME  CUTTING 

X)  CONTINUE 

SIGEST  =  SIGEFF 
SIGOLD  *  SIGEST 
BOLD  =BETAEFF 
ZI  =ZI0 
ZIEST  =  ZIO 
ZD  aZDO 
ALPHA  =  ALPHAO 
AVGEPS  =  AVGEPSO 

RESTORE  OLD  MATERIAL  CONSTANTS 

G  =  GO 
AK  =  AKO 
ZO  «  ZOO 
Z2  =  Z20 
Z3  =  Z30 
A1  =  AlO 
A2  =  A20 
AM2  *  AM20 

UPDATE  NEW  RATE  MATERIAL  CONSTANTS  WITH  SUBINCREMENT  ISUB 

TFACrOR  =  1./  ISUB 

DTSUB  =  DTAU  •  TFACTOR 
DGSUB  =DG  ♦TFACTOR 
DKSUB  =DK  *TFACTOR 
DANSUB  =  DAN  *  TF ACTOR 
DZOSUB  =  DZO  *  TF  ACTOR 
DZ2SUB  =  DZ2  •  TF  ACTOR 
DZ3SUB  =DZ3  'TF ACTOR 
DAISUB  =  DAI  •  TF  ACTOR 
DA2SUB  =DA2  *TFACTOR 
DAM2SUB  =  DAM2  »  TFACTOR 

COMPUTE  DEVITORIC  STRESS  AND  SUBINCREMENTAL  DEVITORIC  AND 
VOLUMETRIC  STRAIN  RATES 

SIGHYD  =  (STRESOd)  +  STRES0(2)  +  STRES0(4))  /  3.0 
DO  20, 1  =  1,4 
C 

EIN(I)  =  EINO(I) 

BETA(I)  =BETA0(I) 

BETAESTd)  =  BETAOd) 
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IF(1^Q.3)THEN 

FACn=0. 

ELSE 

FACT!  =  1. 

ENDIF 

C 

DEVSIG(I)  =  STRESOd)  -  FACTTSIGHYD 
DDEPSd)  =  DEPS(I)  *TFACTOR 
EPS(I)  =  EPSOd) 

C 

20  CONTINUE 


D0200JSUB=1,ISUB 

UPDATE  ALL  TEMPERATURE  DEPENDENT  MATERIAL  CONSTANTS  TO 
END  OF  SUBTIME  INCREMENT  STEP 

G  =  G  +  DGSUB 
AK  =  AK  +  DKSUB 
AN  =  AN  +  DANSUB 
ZO  =  ZO  +  DZOSUB 
Z2  =  Z2  +  DZ2SUB 
Z3  =  Z3  +  DZ3SUB 
A1  =  A1  +  DAISUB 
A2  =  A2  +  DA2SUB 
AM2  =  AM2  +  DAM2SUB 
C 

AVGEPS  =  (EPS(1)+DDEPS(1)  + 

&  EPS(2)+DDEPS<2)  + 

&  EPS(4)+DDEPS(4))  /  3.0 
C 

DO  30 1  =  1,4 

EPSd)  =  EPSd)  +  DDEPSd) 

C 

IFd  .EQ.3)THEN 
DEVEPSd)  =  EPSd) 

ELSE 

DEVEPSd)  =  EPSd)  -  AVGEPS 
ENDIF 

30  CONTINUE 
C 

ICOUNT  =  0 
C 

300  CONTINUE 
C 

ZTOT  =  7IEST  +  ZD 
C 

IF(SIGEST  .LE.  1.E-12)  THEN 
SIGEST=  l.E-12 
DEIEFF  =  0.0 
ELSE 
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XTMPl  s  (ZTar/SIGEST)**2 
XTMP2  » -05*XTMPr*AN 
DEIEFF  =  IXPEXP(XTMP2) 

ENDIF 

C 

XLAM  =  SQRT3  •  DEIEFF/SIGEST 
C 

DEIEFF  =  DEIEFF  *  DTSUB 
C 

SIGNEW  =  0. 

SSUM  *0. 

PWDOT»0. 

C 

DO  40 1=1,4 
IF(I.EQ.3)THEN 
FACTl  =  2 
FACr2  =  0 
FACT3  =  05 
ELSE 

FACTl  =  1 
FACT2  =  1 
FACn  =  l 
ENDIF 

ESTIMATE  PLASTIC  STRAINS  AND  STRESSES 
(ENGINEERING  PLASTIC  SHEAR  STRAINS  ARE  COMPITTED) 

DEPSI(I)  =  XLAM  *  DEVSIG(l)  *  FACTl 

COMPUTE  THERMAL  DIFFERENTIAL  TERMS 

THETAS  =  0. 

C 

EINEST(I)  =  EIN(I)  +  (DEPSKI)  ♦  DTSUB)  +  THETAS 
C 

DEVSIG(I)  =  2.*G*(DEVEPS(I)-EINEST(I))  *  FACTS 
C 

STRESSd)  =  DEVSlCd)  +  FACT2  •  AK  •  AVGEPS 
C 

PWDOT  =  PWDOT  +  (STRESS(I)*DEPSI(I)) 

C 

SIGNEW  =  SIGNEW  +  FACTl  •  DEVSIG(I)«2 
SSUM  =SSUM  +  FACT!  •  STRESS(I)”2 
C 

40  C(>JT1NUE 
C 

SIGNEW  =  SQRT(1.5*S1GNEW) 

SSUM  =SQRT(SSUM) 

C 

IBCOUNT  =  0 

41  BNEW  =  0.0 
C 
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DO  42 1=1,4 
IF(I.EQ.3)THEN 
FACTl  =  2 
ELSE 
FACri  =  1 
ENDIF 
C 

42  BNEW  =BNEW  +  FACTl  *  BETAEST(I)**2 
BNEW  =SQRT(BNEW) 

C 

DO  43 1=1,4 
C 

C  COMPUTE  DRAG  STRESS  VECTORS 
C 

Va)  =  BETAESTd) 

IF(BNEW  .GT.  l.E-15)  V<I)  =  BETAESTd)/  BNEW 
Ud)  =  STRESSd) 

IF(SSUM  .GT.  l.E-15)  Ud)  =  STRESSd)  /  SSUM 
C 

43  CONTINUE 
C 

BSSUMl  =  0. 

BSSUM2  =  0. 

ZD  =0. 

C 

rv^  4A  1-14 

1¥(IEQ3)THEN 

FACT=2 

ELSE 

FACT=1 

ENDIF 

C 

BTERMl  =  (Z3  *  U(I)  -  BETAESTd))  *  PWDOT 
BTERM2  =  Z1  •  Vd)  *((BNEW/Z1)**R2) 

C 

THETA2=  BETAESTd) ’DZ3SUB/Z3 
C 

BETADOTd)  =  AM2’BTERM1  -  A2*BTERM2 
C 

BETAESTd)  =  BET  Ad)  +  BETADOTd)*DTSUB  +  THETA2 
C 

ZD  =  ZD  +  BETAESTd)  »  Ud)  *  FACT 
CBIAXI  BSSUMl  =  BSSUMl  +  BETADOTd)  *  BETADOTd)  *  FACT 
CBIAXI  BSSUM2  =  BSSUM2  +  BETAESTd)  *  BETADOTd)  *  FACT 

44  CONTINUE 
C 

C  SPECIAL  TRAP  FOR  BETA  SNAP  REVERSALS  AT  HIGH  STRAIN  RATES 
C 

ZTOT  =  ZD  +  ZIEST 
IF(ZTOT.LT.O)THEN 
BSCALE  =  ZIEST  /  ABS(ZD) 

ZD  =  ZD*BSCALE*0.5 
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DO  45 1  =  1,4 

45  BETAESTd)  =  BETAEST(I)  *  BSCALE  *  05 
ENDIF 
C 

BERROR  =  DABS{BNEW-BOLD2) 

BMAX  =DMAXl(BNEW,BOLD2) 

IF{BMAX.GT.  TOLER) 

&  BERROR  =  BERROR  /  BMAX 
C 

BOLD2  =  BNEW 
C 

IF  (BERROR  .GT.  TOLER)  THEN 
IBCOUNT  » IBCOUNT  +  1 
IF  (IBCOUNT  -GT.  5)  THEN 
ISUB  =  ISUB*2 
IFdSUB  .GT.  64)  IDBUG  =  1 
IF(ISUB.GT.  128)  THEN 
WRITE  (6,*) '  BETA  REFUSES  TO  CONVERGE ' 

STOP 

ELSE 

GOTO  100 
ENDIF 
ENDIF 
GO  TO  41 
ENDIF 
C 

CBIAXI  IF(BSSUM1  .LT.  1  .OE-12)  BSSUMl  » l.E-12 

CBIAXI  VVDOT  =  BSSUM2  /  BSSUMl 
C 

CBIAXI  IF(ABS(VVDOT).LE.l.)  THEN 
CBIAXI  THETA  =  ACOS(VVDOT) 

CBIAXI  ELSEIF(  VVDOT  .GT.  1 .)  THEN 

CBIAXI  THETA  =  0. 

CBIAXI  ELSE 

CBIAXI  THETA  =  3.141592654 

CBIAXI  ENDIF 
C 

CBIAXI  ALDOT  =  AM2  *  ( ALPH A1  *S1N(THETA)- ALPHA)  *  PWDOT 
C 

CBIAXI  ALEST  =  ALPHA  +  ALDOT  *  DTSUB 
C 

ZTERMl  =  Z1*(ABS(ZIEST-Z2)/Z1)»*R1 
C 

THETAl  =  DZ2SUB 
C 

CBIAXI  ZIDOT  =  AMr(Zl+ALEST*Z3-ZIEST)*PWDOT 
CBIAXI  &  -  A1*ZTERM1 

C 

ZIDOT  =  AMr(Zl-ZIEST)»PWDOT 
&  -  Al’ZTERMl 

C 

C  COMPUTE  UPDATED  ZTOT  VALUE 


53 
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ZIEST  *  ZI  +  ZIDOPDTSUB  +  THETAl 

NOW  INVESTIGATE  CONVERGENCE 

SERROR  =  DABS(SIGEST-SIGNEW) 

SIGMAX=  DMAXKSIGNEW^IGEST) 

IF{  SIGMAX  .GT.  TOLER) 

&  SERROR  =  SERROR  /  SIGMAX 

DEPERR  =  DABS(DEIEFF-DOLD) 

DEPMAX  =»  DMAXl(DEIEFF,DOLD) 

IF(  DEPMAX  .GT.  TOLER) 

&  DEPERR  =  DEPERR  /  DEPMAX 

ZERROR  =  DABS{ZlEST-ZOLD) 

ZMAXs  DMAXKZIEST^OLD) 

IFIZMAX.GT.  TOLER) 

&  ZERROR  =  ZERROR  /  ZMAX 

ERROR  SERROR  +  DEPERR  +  BERROR  +  ZERROR 

IF(IDBUG.EQ.1)THEN 

WRITE(6  *)  •  TIME  KTR ICOUNT  AND  ISUB  ^TIME,  KTR,ICOUNT,ISUB 
WRITE(6/) 

WRITE(6/)  'SIG  EST  NEW  SERR SIGEST^IGNEW^ERROR 
WRITE(6/)  'DEP  OLD  NEW  EPERR  DOLD,DEIEFF,DEPERR 
WRITE(6,*)  ’B  OLD  NEW  BERR BOLD,BNEW,BERROR 
WRITE(6,*)  'Z  OLD  NEW  ZERR  *,  ZOLD^IEST^RROR 
WRrrE(6/) 

ENDIF 

CHECK  FOR  CONVERGENCE 

IF(  ERROR.GT.  TOLER)  THEN 

NONCONVERGENT  SOLUTION  ARRIVES  HERE 

SIGEST  =  SlGOLD*(l  .-RELAX)  +  SIGNS W  *  (RELAX) 

SIGOLD  =  SIGEST 
DOLD  =  DEIEFF 
BOLD  =  BNEW 
ZOLD  =  ZIEST 

UPDATE  AND  CHECK  CONVERGENCE  COUNT 
ICOUNT  =  ICOUNT+1 
IF  (ICOUNT  .GT.  30)  THEN 
MAKE  ADDITIONAL  TIME-STEP  CUTS 
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ISUB  =  ISUB  •  2 
C 

IF  {ISUB  .EQ.  128)  IDBUG  *  1 
C 

C  TERMINATE  NONCONVERGE  SOLUTIONS 
C 

IF  (ISUB  .GT.  128)  THEN 

WRITE(6,*)  •  SOLUTION  REFUSED  TO  CONVERGE ' 

STOP 

ENDIF 

C 

C  RESTART  ITERATION  SCHEME  WITH  NEW  TIME-STEP  CUT 
C 

GOTO  100 
C 

ENDIF 

C 

GOTO  300 
C 

ELSE 

C 

C  CONVERGENCE  SOLUTION  ARRIVES  HERE 
C 

ALPHA  =  ALEST 
ZI  =ZIEST 
C 

DO  501  =  1,4 
EIN(I)  =EINEST(I) 

BETA(I)  =BETAEST(I) 

50  CONTINUE 
C 

ENDIF 

C 

C  END  OF  SUBTIME  LOOPS 
C 

200  CONTINUE 
C 

SIGEFF  =  SIGEST 
BETAEFF  =  BNEW 
C 
C 

C  DEBUG  OUTPUT  STATEMENTS 
C 

CJOE2  IF  (IPT  .EQ.  1  .AND.  KTR  .EQ.  INTER)  THEN 
C 

CJOE2  WR1TE(6,*) '  TIME  KTR  ICOUNT  AND  ISUB  ^TIME,  KTR.lCOUNT,ISUB 
CIOE2  WRITE(6,*) 

CJOE2  WRITE(6,*)  'SIG  EST  NEW  SERROR  SIGEST^IGNEW^ERROR 
CJOE2  WR1TE(6,*)  'DEP  OLD  NEW  EPERROR  DOLD,DEIEFF,DEPERR 
CJOE2  WR1TE(6,*)  B  OLD  NEW  BERROR  *,  BOLD,BNEW,BERROR 
CJOE2  WRITE(6,*)  Z  OLD  NEW  ZERROR Z0LD;ZIEST,ZERR0R 
CJOE2  WRITE(6,*) 
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CJOE2  WRITE(6  *) '  FOR  ELEMENT  AND  G  POINT NEL,  IPT 
CJOE2  WRrrE(6  •) 

CJOE2  WR1TE(6  *) '  EINEST  EINEST 

CJOE2  WR1TE(6,*)  ’  STRESS  STRESS 

CJOE2  WRITE(6  *) '  BETAEST  *,  BETAEST 

CJOE2  WRITE(6,*) '  ZIEST,  DEIEFF,  ZD ZIEST,  DEIEFF,  ZD 

C 

CJOE2  ENDIF 
C 

RETURN 

END 

C 

SUBROUTINE  STGET2(ARRAY) 

C 

IMPLICIT  DOUBLE  PRECISION  ( A-H.O*Z ) 

C 

REAL’S  ARRAY(60) 

C 

COMMON  /BPSTAT2/  E1N(4),  ZI,S1GEFF,  BETAEFF,  EPEFF,  ALPHA, 
&  ZD,  BETA(4),  EPS0(4) 

C 

EIN(l)  =ARRAY(1> 

EIN(2)  =  ARRAY(2) 

EIN(3)  =  ARRAYO) 

EIN(4)  =  ARRAY(4) 

BETA(l)  =  ARRAY(5) 

BETA(2)  =  ARRAY(6) 

BETA(3)  =  ARRAY<7) 

BETA(4)  =  ARRAY(8) 

ZI  =  ARRAY(9) 

SIGEFF  =  ARRAY(IO) 

EPEFF  =  ARRAYO  1) 

BETAEFF=  ARRAYO  2) 

ALPHA  =ARRAY03) 

ZD  =  ARRAY04) 

EPS0O)  =  ARRAYO5) 

EPS(X2)  =  ARRAY06) 

EPS0(3)  =  ARRAYO?) 

EPS0(4)  =  ARRAYO8) 

C 

RETURN 

END 

C 

SUBROUTINE  STPUT2(ARRAY) 

C 

IMPLICIT  DOUBLE  PRECISION  (  A-H,0-Z  > 

C 

REAL’S  ARRAY(60) 

C 

COMMON  /BPSTAT2/  EIN(4),  ZI,  SIGEFF,  BETAEFF,  EPEFF,  ALPHA, 
&  ZD,BETA(4),EPS0(4> 

C 
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ARRAY(1)»EIN(1) 
ARRAY(2)«EIN(2) 
ARRAYO)  =  EIN(3) 
ARRAY(4)  a  EIN(4) 
ARRAY(5)  a  BETA(l) 
ARRAY(6)  a  BETA(2) 
ARRAY(7)aBETA(3) 
ARRAY(8>aBETA(4) 
ARRAY(9>aZI 
ARRAY(10)aSlGEFF 
ARRAY(11)»  EPEFF 
ARRAY(12)aBETAEFF 
ARRAY(13)a  ALPHA 
ARRAY{t4)aZD 
ARRAY(15)a  EPS0(1) 
ARRAY(16)aEPS0(2) 
ARRAY(17)aEPS0(3) 
ARRAY(18)a  EPS0{4) 

C 

RETURN 

END 

C 
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Appendix  B 


User-Defined  Subroutines  for  Three-Dimensional  Brick  Elements 

(file  cuser3_dbeta.f) 


SUBROUTINE  CUSER3  (NG>IEL,IPT,STRESS,EPS,STRAIN,DEPS,DEPST, 

A  THSTR1,THSTR2,KTR,INTEROT.ARRAY,IARRAY,D, 

B  ALFA,CTD,ALFAA,CTDD,CTI,TMP1,TMP2,T1ME,DTAU, 

C  PHIST,PRST,PH1STUDPSP,INTEG,ISUBM,INDNL,KEY) 

C*I 

ei 

C*I 

C*I  THIS  SUBROUTINE  IS  TO  BE  SUPPLIED  BY  THE  USER  TO  CALCULATE 
C*I  THE  STRESSES  AND  CONSTITUTIVE  MATRIX  OF  A  SPECIAL  MATERIAL 
C*I 

C*I  THIS  SUBROUTINE  IS  CALLED  IN  USER2  FOR  EACH  INTEGRATION  POINT 
C*I  FOR  EACH  3-D  SOLID  ELEMENT  TO  PERFORM  THE  FOLLOWING  OPERATIONS : 

C*I 

C*I  KEY.EQ.1  INITIALIZE  THE  WORKING  ARRAYS  DURING  INPUT  PHASE 
C*I 

ei  KEY.EQ.2  CALCULATE  ELEMENT  STRESSES 
C*I 

C*I  KEY.EQ.3  CALCULATE  THE  STRESS/STRAIN  MATRIX 
C*! 

C*I  KEY.EQ.4  PRINT  CALCULATED  STRESSES  AND  OTHER  DESIRED 
ei  VARIABLES  DURING  STRESS  PRINT-OUT 

ei 

C*I 

C*I  THE  FOLLOWING  VARIABLES  ARE  USED  TO  PERFORM  THE  ABOVE  OPERATIONS : 
C*I 

C*!  NG  ELEMENT  GROUP  NUMBER 

OI 

C*I  NEL  ELEMENT  NUMBER 

C*I 

C*I  IPT  INTEGRATION  POINT  NUMBER 

ei 

C*I  STRESS(6)  STRESS  COMPONENTS,  RECEIVED  AT  TIME  TAU 
C'l  AND  UPDATED  BY  USER-SUPPLIED  CODING  TO 

CM  CORRESPOND  TO  TIME  TAU+DTAU 

CM 

CM  EPS(6)  TOTAL  STRAIN  COMPONENTS  AT  TIMET. 

C*I  IN  CASE  OF  LARGE  STRAIN  FORMULATION 

CM  (U.L.H.),  EPS{4)  ARE  ELASTIC  STRAINS  AT 

CM  TIME  T. 

CM 

CM  STRAIN(6)  TOTAL  STRAIN  COMPONENTS  AT  TIME  T+DT. 

CM  INCASEOFU.L.H.: 

CM  A)  FOR  KEY.EQ.2  -  ELASTIC  STRAINS 

C*I  AT  TIME  TAU+DTAU  (FOR  KTR.EQ.O  - 
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C*I 

ei 

oi 

c*i 

ei 

c*i 

c*i 

c*i 

CT 

c*i 

c*i 

c*i 

ei 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

c*i 

C*I 

C*I 

C*I 

ei 

c*i 

c*i 

c*i 

ei 

ei 

ei 

c*i 

c*i 

ei 

c»i 

c*i 

c*i 

c*i 

ei 

c*i 

C‘I 

cn 

c*i 

c*i 

c*i 

c*i 

ei 

c*i 

c*i 


TRIAL  ELASTIC  STRAINS) 

B)  FOR  KEY.GT.2  -  ELASTIC  STRAINS  AT 
TIME  T+DT 

DEPS(6)  SUBDIV IDED  INCREMENTAL  STRAIN  COMPONENTS 
DEPS(1)  =  ( STRAlN(l)  -  EPS(I)  )/INTER 
-  DEPSTd) 

( PASSED  TO  SUBROUTINE  CUSER2  BY  THE 
PROGRAM  ADINA ) 

DPSP(6>  INCREMENT  OF  INELASTIC  STRAIN  (PLASTIC 
AND/OR  CREEP  AND/OR  VISCOPLASTIC,  ETC.) 

IN  THE  SUBDIVISION.  CALCULATED  BY  USER- 
SUPPLIED  CODING  AND  PASSED  TO  THE  PROGRAM 
ADINA.  USED  FOR  U.L.H.  ONLY. 

DEPST<6)  COMPONENTS  OF  SUBINCREMENTAL  THERMAL 
STRAIN 

( PASSED  TO  SUBROUTINE  CUSER2  BY  THE 
PROGRAM  ADINA ) 

THSTR1(6)  TOTAL  THERMAL  STRAIN  AT  TIME  TAU 

FOR  KEY.EQ.4  -  THERMAL  STRAIN  AT  TIME  T 

THSTR2(6)  TOTAL  THERMAL  STRAIN  AT  TIME  TAU+DTAU 

FOR  KEY  EQ.4  -  THERMAL  STRAIN  AT  TIME  T+DT 

INTER  NUMBER  OF  SUBDIVISIONS  FOR  THE  STRAIN 

INCREMENTS  (INTER  =  INT(PROP(123)) 

KTR  CURRENT  SUBDIVISION  NUMBER 

EQ.0  CALCULATION  OF  TRIAL  ELASTIC 
STATE,  IN  CASE  INTEG=1 
EQ.l  FOR  HRST  SUBDIVISION 
EQ.INTER  FOR  LAST  SUBDIVISION 

SCP(4)  SOLUTION  CONTROL  PARAMETERS 

ARRAY(60)  WORKING  ARRAY  (REAL) ,  RECEIVED  AT 

TIME  TAU  AND  UPDATED  BY  USER-SUPPLIED 
CODING  TO  CORRESPOND  TO  TIME  TAU+DTAU 

IARRAY(2)  WORKING  ARRAY  (INTEGER) ,  RECEIVED  AT 

TIME  TAU  AND  UPDATED  BY  USER-SUPPLIED 
CODING  TO  CORRESPOND  TO  TIME  TAU+DTAU 

D(6.6)  STRESS/STRAIN  MATRIX ,  TO  BE  CALCULATED 

BY  USER-SUPPLIED  CODING 

ALFA  COEFFICIENT  OF  THERMAL  EXPANSION  AT 

TIME  TAU 


CTD(5)  TEMPERATURE-DEPENDENT  MATERIAL  CONSTANTS 

AT  TIME  TAU 
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C*I 

ei 

ALFAA 

OOEFFIOIENT  OF  THERMAL  EXPANSION  AT 

ei 

TIME  TAU+DTAU 

C*I 

C*I 

OTDEKS) 

TEMPERATURE-DEPENDENT  MATERIAL  OONSTAN’R 

C*I 

AT  TIME  TAU+DTAU 

ei 

C*I 

on(8) 

TEMPERATURE-INDEPENDENT  MATERIAL  OONSTANTS 

C*! 

C*I 

TMPl 

TEMPERATURE  AT  INTEGRATION  POINT  IPT  AT 

C*I 

TIME  TAU.  FOR  KEY.EQ.4  -  TEMPERATURE  AT 

ei 

TIMET 

C*I 

ei 

TMP2 

TEMPERATURE  AT  INTEGRATION  POINT  IPT  AT 

C*I 

TIME  TAU+DTAU.  FOR  KEY.EQ.4  -  TEMPERATURE 

0*1 

ATTIMET+DT 

0*1 

ei 

TIME 

TIME  AT  OURRENT  STEP ,  T+DT 

0*1 

0*1 

DT 

TIME  STEP  INOREMENT ,  D T 

ei 

0*1 

PHIST(3,3) 

MATRIX  OONTAINING  DIREOTION  OOSINES  OF 

0*1 

PRINOIPAL  STRETOH  DIRECTIONS,  IN  CASE  OF 

0*1 

U.L.H.  FORMULATION 

0*1 

0*1 

PRST(3) 

PRINCIPAL  STRETCHES,  IN  CASE  OF  U.L.H. 

0*1 

0*1 

PHIST1(3) 

DIRECTION  COSINES  OF  THE  FIRST  PRINCIPAL 

0*1 

STRETCH  (STRETCH  WHICH  HAS  MAXIMUM 

0*1 

DEVIATION  FROM  1.0) 

0*1 

0*1 

INTEG 

INTEGRATION  PARAMETER  FOR  STRESS 

0*1 

INTEGRATION 

ei 

EQ.O  -  FORWARD  INTEGRATION 

0*1 

EQ.l  -  BACKWARD  INTEGRATION 

0*1 

(RETURN  MAPPING) 

0*1 

0*1 

ISUBM 

FLAG  FOR  CONTINUATION  OF  SUBDIVISION 

0*1 

IN  TIME  STEP,  APPLICABLE  FOR  INTEG.EQ.l 

0*1 

EQ.O  -  CONTINUATION 

0*1 

EQ.-l  -  STOP  OF  SUBDIVISION 

0*1 

IN  THE  USER-SUPPLIED  CODING  THE  FLAG 

0*1 

(INITIALLY  EQ.O)  MUST  BE  SET  TO  -1 

0*1 

WHEN  CRITERIA  FOR  STOPPING  SUBDIVISIONS 

0*1 

ARE  REACHED 

ei 

0*1 

INDNL 

FLAG  FOR  ELEMENT  FORMULATION  (NPAR(3)) 

0*1 

EQ.l  -  MATERIALLY  NONLINEAR  ONLY  (M.N.O.) 

0*1 

EQ.2  -  LARGE  DISPLACEMENTS  AND  SMALL 

0*1 

STRAINS  (T.L.) 

0*1 

EQ.3  -  LARGE  DISPLACEMENTS  AND  LARGE 

0*1 

STRAINS  (U.L.H.) 

0*1 

0*1 
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C*I  NOTE  THAT  THE  VARIABLES  PASSED  TO  THE  SUBROUTINE  WHEN  KEY=3,4 
C*1  ARE  THESE  CALCULATED  LAST:  I.E.,  CALCULATED  IN  THE  LAST 
C*I  SUBDIVISION  WHEN  KEY=2.  HENCE  THESE  VARIABLES  ARE  NOT 
C*I  CALCULATED  WHEN  KEY=3,4. 

C*I 

ei 

C*l 

IMPLICIT  DOUBLE  PRECISION  ( A-H.O-Z ) 

C 

DIMENSION  STRESS(6)OTAIN(6),DEPS(6),ARRAY(60),IARRAY(2),D(6^), 

A  CTD(5),CTI(8),EPS(6),SCP(4) 

DIMENSION  CTDD(5),DEPST(6)JHSTR1  (6).THSTR2(6) 

DIMENSION  PHIST(3^),PRST(3XPHIST1  (3),DPSP(6) 

C 

C 

REAL*8  DTAU 

C  TIME  INCREMENT  STEP,  DT 

REAL*8  E 

C  ELASTIC  MODULUS 

REAL*8  G 

C  SHEAR  MODULOUS 

INTEGER  HINT 

C  INTEGRATION  POINT  WRITTEN  TO  HISTORY  FILE 

INTEGER  HELE 

C  ELEMENT  NUMBER  ASSOCIATED  WITH  HIPT 

INTEGER  I 

C  DO  LOOP  COUNTER 

INTEGER  J 

C  DO  LOOP  COUNTER 

C 

C  DIRECTIONAL  BODNER-PARTOM  MATERIAL  CONSTANTS: 

C 

C  TEMPERATURE  INDEPENDENT  CONSTANTS 
C 

REAL*8  AMI 
REAL’S  Z1 
REAL’S  R1 
REAL’S  R2 
REAL’S  DO 
REAL’S  ALPHAl 
C 

C  TEMPERATURE  INDEPENDENT  CONSTANTS 
C 

REAL’S  AN 
REAL’S  AM2 
REAL’S  ZO 
REAL’S  Z2 
REAL’S  Z3 
REAL’S  A1 
REAL’S  A2 
C 

C  TEMPERATURE  DIFFERENTIALS 
C 

REAL’S  DAM2 
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REAL‘8  DAN 
REAL’S  DZO 
REAL’S  DZ2 
REAL’S  DZ3 
REAL’S  DAI 
REAL’S  DA2 
C 

COMMON /MATCONST/  E,  C,  AK.  DC,  DK, 

&  AMI,  AM2,  ALPHAl,  DAM2. 

&  Zl,  Z3,  Rl,  R2,  DO, 

&  AN,  ZO,  Z2,  Al,  A2, 

&  DAN,  DZO,  DZ2,  DZ3,  DAI,  DA2 

C 

C  COMMON  BLOCK  FOR  STATE  VARIBLES 
C 

COMMON  /BPSTAT3/  EIN(6),  ZI,  SICEFF,  BETAEFF,  EPEFF,  ALPHA, 

&  ZD,  BETA(6),  EPS0(6) 

C 

C  RECALL  STATE  VARIBLES  AT  TIME  TAU 
C 

CALL  STGET3{ARRAY) 

C 

GOTO  (1,2,3,4),  KEY 
C*! 

C’l 

C’l  KEY  =  1 
C*l 

ei  INITIALIZE  COMPONENTS  OF  REAL  AND  INTEGER  WORKING  ARRAYS 
C’l  ( INITIALIZE  ARRAY(60)  AND  IARRAY(2) ) 

C’l 

1  CONTINUE 
C’l 

C’l  ’”  INSERT  USER-SUPPLIED  CODING 

C’l 

C’l 

C  INITIALIZE  STATE  VARIBLES 
C 

EIN(l)  =0.0 
EIN(2)  =  0.0 
EIN(3)  =0.0 
E1N(4)  =0.0 
EIN(5)  =  0.0 
E1N(6)  =  0.0 
BETA(1)=0.0 
BETA(2)  =  0.0 
BETA(3)  =  0.0 
BETA(4)  =  0.0 
BETA(5)  =  0.0 
BETA(6)  =  0.0 
ZI  =  CTD{3) 

SIGEFF  =  0.0 
EPEFF  =0.0 
BETAEFF  =  0.0 
ALPHA  =0.0 
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ZD  =  0.0 


C 

C  STORE  STATE  VARIBLES 
C 

CALL  STPUT3(ARRAY) 

C 

RETURN 

C*1 

C*I 

C*I  KEY  =  2 
C*I 

C*I  INTEGRATION  OF  ELEMENT  STRESSES 
C*I  { CALCULATE  STRESS(6)) 

C*I 

2  CONTINUE 
C‘I 

C*I  INSERT  USER-SUPPLIED  CODING 

ei 

ei 

C  BODNER-PARTOM  MATERIAL 
C 

C  FOR  3-D  SOLID  ELEMENT 
C 

C  GET  CURRENT  VALUES  OF  MATERIAL  CONSTANTS  FOR  GIVEN  TEMPERATURE 

C 

C 

AMI  =CTI(1) 

Z1  =  CTI(2) 

R1  =  CTI(3) 

R2  =  CTI(3) 

DO  =  CTI(4) 

ANU  =CTI(5) 

ALPHAl  =  CTI{8) 

C 

A1  =  CTI(6)  *  EXP(-CTI(7)/(TMP1  +273.)) 

A2  =A1 

DAI  =  CTI(6)  *  EXP(-CTI(7)/(TMP2+273.))-Al 
DA2  =  DAI 
C 

G  =  CTD(1)/(2.*(1.+ ANU)) 

DG  =  CTDD(])/(2.»(1.+  ANU))  -  G 
C 

AK  =  CTD(l)/(I.-2.*  ANU) 

DK  =  CTDD(l)/(l.-2.*  ANU)  -  AK 
C 

AN  =  CTD(2) 

DAN  =  CTDD(2)-CTD(2) 

C 

Z2  =  CTD(3) 

DZO  =  CTDD(3)-CTD(3) 

DZ2  =  CTDD(3)-CrD{3) 

C 

Z3  =  CTD(4) 

DZ3  =  CTDD(4)-CTD(4) 
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c 

AM2=  CTD(5) 

DAM2  =  CTDD{5)-CTD(5) 

C 

RELAX  =  SCP(1) 

TOLER  =»SCP{2) 

HINT  « inX(SCP(3)) 

HELE  =  IF1X(SCP(4)) 

C 

C  COMPUTE  INCREMENTAL  STRAINS  WITH  INCREMENTAL  THERMAL  STRAINS 
C 

IF(KTR.EQ.1)THEN 
00201  =  1,6 

20  EPSOIl)  =  EPS(I)  -  THSTRl  (I) 

ENDIF 

C 

DTAUl  =  DTAU  /  FLOAT!  INTER  ) 

C 

C  CALCULATE  INCREMENTAL  STRESSES  AND  UPDATE  STRESS  VECTOR 
C 

CALL  DBODNER3  ( DEPS,  STRESS,  DTAUl, 

+  NEL,  IPT,  KTR,  TIME,  INTER,  RELAX,  TOLER ) 

C 

DO  30  I  =  1,6 

30  EPS0(I)  =  EPS0{I)  +  DEPS(l) 

C 

c - 

C  STORE  STATE  VARIBLES  AT  TIME  TAU  + DTAU 
C 

CALL  STPUT3(ARRAY) 

C 

RETURN 

ei 

c*i 

C*I 

C»I  KEY  =  3 
C*I 

C*I  FORM  CONSTITUTIVE  LAW 
C*I  ( CALCULATE  D(6,6) ) 

C*I 

3  CONTINUE 
C*I 

C*I  INSERT  USER-SUPPLIED  CODING 
C*I 
C 
C 

C  THE  TEMPERATURE  DEPENDENT  ELASTIC  MATRIX  IS  DETERMINED 
C  WITHOUT  A  PLASTIC  STRAIN  CORRECTION 
C 

E  =CTDD(1) 

ANU  =  CTI(5) 

C 

DDT  =  DTAU/FLOAT(INTER) 

C 
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DO  315  I  =  U 
D0315J  =  1^ 

315  CKU)  =  0.D0 
CM  =  E/(1.-2.*ANU) 

AE  =  (1.  + ANU)/E 
CP  =  AE 
C 

C  IN  CASE  CORRECTION  FOR  VISCOPLASTIC  FLOW  IN  THE  TIME  STEP, 
C  CORRECT  THE  CONSTANT  CP 
C 

CJOE  IF  (EST.GT.YLD)  CP  =  CP  +  1.5*DDT*GAMA*(EST/YLD  -  l.)/EST 
CP  =  1./CP 
Cll=(CM  +  2.*CP)/3. 

C12  =  (CM.CP)/3. 

D(U)  =  C11 
0(1,2)  =  02 
D<l,3)  =  a2 
D(2,2)  =  C11 
D(2.3)  =  02 
D(3,3)  =  01 
D(4.4)  =  0.5*CP 
D(5,5)  =  D(4,4) 

D(6,6)  =  D(5,5) 

003201  =  U 
DO320J  =  I,3 
320  D(J,I)  =  D(U) 

C 

RETURN 

C*l 

C*l 

C*I  KEY  =  4 
OI 

C*I  PRINTING  OF  ELEMENT  RESPONSE 
CM  ( PRINT  STRESS(6),STRAIN(6) ) 

CM 

4  CONTINUE 
CM 

CM  INSERT  USER-SUPPLIED  CODING 

CM 

CM 

IF((NEL.EQ.HELE).AND.(IPT.EQ.HINT))THEN 
7006  FORMAT(2(lX,F8.1),12{lX,ElL4)) 

ZMECHS  =  STRAIN(l)  -  THSTR2(1) 

ZTOT  =ZI+ZD 

WRITE(53,7006)  TIME,TMP2,STRESS(1  ),STRESS{2), 

+  stress(3)z:tot,zi,zd,zmechs,sigeff 

ENDIF 

C 

C  PRINT  HEADING  AND  ELEMENT  NUMBER 
C 

IF  (NEL.EQ.l  .AND.  IPT.EQ.l)  THEN 
WRITE  (6,9000) 

ENDIF 

IF  (IPT.EQ.1)  WRITE  (6,9002)  NEL 
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c 

C  PRINT  STRESSES 
C 

WRITE  (6,9003)  IPT3TRESS,SIGEFF.STRAIN,TMP2 
C 

C  FORMAT  STATEMENTS 
C 

9000  FORMAT  (//.4X3HNEL,2X,3HIPT,6X,9HSTRESS-XX,6X, 

1  9HSTRESS-YY,6X,9HSTRESS*ZZ,6X,9HSTRESS-XY,6X, 

2  9HSTRESS-X2,6X,9HSTRESS-YZ,4X, 

2  IIHEFF  STRESS,/ 

3  17X,9HSTRAIN-XX,6X, 

1  9HSTRAIN-YY,6X,9HSTRA1N-ZZ.6X,9HSTRA1N-XY,6X, 

2  9HSTRAIN-XZ,6X,9HSTRA1N-YZ,4X, 

2  IIHTEMPERATURE,  /) 

9002 FORMAT!/ 17  /) 

9003  FORMAT  (7X,I5,7(2X,E13.6)/13X.7(2X,E13.6),/> 

9004  FORMAT  (7X,I5.8(2X,E13.6)/27X,3(2X,E13.6)30X,2(2X,E13.6)  / 

1  27X,3(2X,E13.6)/) 

C 

CALL  STPUT3(ARRAY) 

C 

RETURN 
ORLE  END 
END 
C 

SUBROUTINE  DBODNER3(DEPS,  STRESS,  DTAU, 

+  NEL,  IPT,  KTR,  TIME,  INTER,  RELAX,  TOLER) 

C 

IMPLIQT  DOUBLE  PRECISION  ( A-H,0-Z ) 

C - 

c 

REAL'S  DEVEPS(6) 

REAL'S  EPS(6) 

C 

REAL'S  AVGEPS 

C  AVERAGE  NORMAL  STRAIN  INCREMENT  AT  TIME  T 

REAL'S  DDEPS(6) 

C  SUBINCREMENTIALDEVITORIC  STRAINS 

REAL'S  DEIEFF 

C  EFFECTIVE  INELASTIC  STRAIN  INCREMENT 

REAL'S  DE(^6) 

C  INCREMENTAL  TOTAL  STRAIN  VECTOR 

REAL'S  DEPSI(6) 

C  *  NELASTIC  PORTION  OF  THE  DEVI  ATORIC  STRAIN 

C  VECTOR 

REAL'S  DEVSIG(6) 

C  DEVIATORIC  STRESS  VECTOR  OF  THE  STRESS  STATE 

C  AT  TIMET 

INTEGER  ICOUNT 

C  ITERATION  LCX)P  COUNTER 

REAL'S  FAcrri 
C  INTEGER  FACTOR 

REAL'S  FACT2 
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C  INTEGER  FACTOR 

INTEGER  INTER 

C  NUMBER  OF  SUBINCREMENTS 

INTEGER  IPT 
C  GAUSS  POINT 

INTEGER  ISUB 

C  NUMBER  OF  SUBINCREMENTS  WITHIN  C 

INTEGER  KTO 

CURRENT  SUBDIVISION  NUMBER 
EQ.  1  FOR  FIRST  SUBDIVISICMsI 
EQ.  2  FOR  LAST  SUBDIVISION 
INTEGER  NEL 

ELEMENT  NUMBER 
REAL*8  RELAX 

RELAXATION  FACTOR  FOR  NEW  STRESS  ESTIMATE 
REAL*8  TOLER 

TOLERANCE  FOR  CONVERGENCE 
REAL'S  SIGNEW 

NEW  EFFECTIVE  STRESS 
REAL'S  SIGOLD 

INITIAL  EFFECTIVE  STRESS 
REAL'S  SIGEST 

EFFECTIVE  STRESS  FOR  TIME  STEP  SIZE  OF  DT/JINT 
REAL'S  SIGHYD 

HYDROSTATIC  STRESS 
REAL'S  SIGMAX 

MAXIMUM  ABSOLUTE  VALUE  OF  EITHER  SIGEFFOR  SIGEFZ 
FOR  INVESTIGATING  CONVERGENCE 
REAL'S  STRES0(6) 

STRESS  VECTOR  AT  TIME  TAU 
REAL'S  STRESS(6) 

STRESS  VECTOR  AT  TIME  TAU+DTAUT 
REAL'S  TIME 

CURRENT  TIME  VALUE 
REAL'S  EINEST(6) 

ESTIMATED  PLASTIC  STRAINS 
REAL'S  EIN(K6) 

PLASTIC  STRAINS  AT  TIME  TAU 
REAL'S  BETADOT(6) 

BETA  RATE  VECTOR 
REAL'S  BETAEST(6) 

ESTIMATED  BETA  VECTOR 
REAL'S  BETA0(6) 

BETA  VECTOR  AT  TIME  TAU 
REAL'S  U(6),  V(6) 

DIRECTIONAL  VECTORS 
REAL'S  ALEST 

ESTIMATED  ALPHA 
REAL'S  ZIEST 

ESTIMATED  ISOTROPIC  DRAG  STRESS 


COMMON  /MATCONST/  E,  G,  AK,  DC,  DK, 
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& 

& 
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AM1,AM2,  ALPHA!,  DAM2, 

Zl,  Z3,  Rl,  R2,  DO, 

AN,  ZO,  Z2,  Al,  A2, 

DAN,  DZO,  DZ2,  DZ3,  DAI,  DA2 
C 

COMMON  /BPSTAT3/  EIN(6),  ZI,  SICEFF,  BETAEFF,  EPEFF,  ALPHA, 
&  ZD,  BETA(6),  EPS0(6) 

C 

SQRT3»SQRT(3.) 

STORE  OLD  STATE  VARIBLES 

DO  10 1=1,6 
STRES0(1)  =  STRESS(I) 

EIN0(I)  =EIN(I) 

BETA0{I)  =  BETA(I) 

I  CONTINUE 

ZIO  =ZI 
ALPHAO  =  ALPHA 
ZDO  =ZD 

STORE  OLD  MATERIAL  CONSTANTS 

GO  «G 
AKO  =AK 
ZOO  =Z0 
Z20  =Z2 
Z30  =Z3 
AlO  =  Al 
A20  =  A2 
AM20  =  AM2 

INITIALIZE  OTHER  VARIBLES 


ISUB  =1 
IDBUG  =  0 

INITIALIZE  VARIABLES  FOR  SUB-TIME  CUTTING 

0  CONTINUE 

SIGEST  =  SIGEFF 
SIGOLD  =  SIGEST 
BOLD  =  BETAEFF 
ZI  =ZI0 
ZIEST  =ZI0 
ZD  =ZD0 
ALPHA  =  ALPHAO 
AVGEPS  =  AVGEPSO 

RESTORE  OLD  MATERIAL  CONSTANTS 

G  =  GO 
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AK  =  AKO 
ZO  =  ZOO 
Z2  =  Z20 
Z3  =  Z30 
A1  =  AlO 
A2  =  A20 
AM2  =  AM20 

UPDATE  NEW  RATE  MATERIAL  CONSTANTS  WITH  SUBINCREMENT  ISUB 

TFACrOR  =  l./ISUB 

DTSUB  =DTAU*TFACTOR 
DGSUB  aDG  *TFACTOR 
DKSUB  =DK  *TFACTOR 
DANSUB  =  DAN  *  TFACTOR 
DZOSUB  =  DZO  *  TFACTOR 
DZ2SUB  =  DZ2  *  TFACTOR 
DZ3SUB  =  DZ3  *  TFACTOR 
DAISUB  =DA1  ’TFACTOR 
DA2SUB  =  DA2  ’  TFACTOR 
DAM2SUB  =  DAM2  *  TFACTOR 

COMPUTE  DEVITORIC  STRESS  AND  SUBINCREMENT AL  DEVITORIC  AND 
VOLUMETRIC  STRAIN  RATES 

SIGHYD  =  (STRESOd)  +  STRES0(2)  +  STRESOO))  /  3.0 
DO  20, 1  *1,6 

EIN(I)  =  EIN0{1) 

BETA(I)  =BETA0(I> 

BETAESTd)  =  BETAOd) 

IF(I.GT.3)THEN 

FACn*0. 

ELSE 

FACTl  *  1. 

ENDIF 

DEVSIG(I)  =  STRESO(I)  -  FACTl  *SIGHYD 
DDEPSd)  =  DEPS(l)  ’TFACTOR 
EPSd)  =  EPSOd) 

0  CONTINUE 


DO  200  JSUB=1,ISUB 

UPDATE  ALL  TEMPERATURE  DEPENDENT  MATERIAL  CONSTANTS  TO 
END  OF  SUBTIME  INCREMENT  STEP 


G  *G  +  DGSUB 
AK  =  AK  +  DKSUB 
AN  =  AN  +  DANSUB 
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20  ®Z0  +  DZOSUB 
22  =  Z2  +  DZ2SUB 
A1  »  A1  +  DAISUB 
A2  =  A2  +  DA2SUB 
AM2  «  AM2  +  DAM2SUB 
C 

AVGEPS  =  (EPS(1)+DDEPS(1)  + 

&  EPS(2)+DDEPS(2)  + 

&  EPS(3)+DDEPS(3))  /  3.0 
C 

DO30I*U 

EPS(I)  =  EPS(I)  +  DDEPS(I) 

C 

IF(1.GT.3)THEN 
DEVEPS(I)  =  EPS(I) 

ELSE 

DEVEPSd)  =  EPS(5)  -  AVGEPS 
ENDIF 

30  CONTINUE 
C 

ICOUNT  =  0 

c 

300  CONTINUE 
C 

ZTOT«ZIEST  +  ZD 
C 

IF(SIGEST  .LE.  lJE-12)  THEN 
SIGEST=1.E-12 
DEIEFF  =  0.0 
ELSE 

XTMPl  =(ZTOT/SlGEST)*^ 
XTMP2  =  -0.5*XTMP1**AN 
DEIEFF*  IXPEXP(XTMP2) 

ENDIF 

C 

XLAM  *  SQRT3  *  DEIEFF/SIGEST 
C 

DEIEFF  «  DEIEFF  *  DTSUB 
C 

SIGNEW  =  0. 

SSUM  =0. 

PWDOT  =  0. 

C 

DO  40 1=1,6 
IF(I.GT.3)THEN 
FACri  =  2 
FACr2=:0 
FACT3  =  1 
ELSE 
FACTl  =  1 
FACT2  =  1 
FACT3  =  2 
ENDIF 
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C  ESTIMATE  PLASTIC  STRAINS  AND  STRESSES 
C  (ENGINEERING  PLASTIC  SHEAR  STRAINS  ARE  COMPUTED) 
C 

DEPSI(l)  =  XLAM  •  DEVSIG(I) '  FACTl 
COMPUTE  THERMAL  DIFFERENTIAL  TERMS 


THETA3  =  0. 


EINESTCI)  =  EIN(I)  +  (DEPSKI)  *  DTSUB)  +  THETAS 

DEVSIG(I)  =  FACn*G*(DEVEPS(I)-  EINESTCI)) 

STRESSCI)  =DEVSIG(I)  +  FACT2*AK'AVGEPS 

PWDOT  =  PWDOT  +  (STRESS(I)*DEPSI(I)) 

SIGNEW  =  SIGNEW  +  FACTl  *  DEVSIG(l)-2 
SSUM  =SSUM  +  FACTl  *  STRESS(I)*»2 

,0  CONTINUE 

SIGNEW  =  SQRT(1.5*SIGNEW) 

SSUM  =SQRT(SSUM) 

IBCOUNT  =  0 
H  BNEW  =  0.0 

DO  42  I=L6 
IF(I.GT.3)THEN 
FACTl  =  2 
ELSE 
FACTl  =  1 
ENDIF 

2  BNEW  =BNEW  +  FACTl  *  BETAEST(I)**2 
BNEW  =SQRT(BNEW) 

DO  43 1*1,6 

COMPUTE  DRAG  STRESS  VECTORS 
V(I)  =  BETAEST(I) 

IFCBNEW  .GT.  l.E-15)  V(I)  =  BETAEST(l)/  BNEW 
U(l)  *  STRESSCI) 

IFCSSUM  .GT.  l.E-15)  U(l)  =  STRESSCI)  /  SSUM 
C 

43  CONTINUE 
C 

BSSUMl  =  0. 

BSSUM2  =  0. 

ZD  =0. 

C 

D044I  =  1,6 
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IF(I.GT.3)THEN 

FACr=2 

ELSE 

FACr=l 

ENDIF 

C 

BTERMl  =  (Z3  *  U(I>  -  BETAESTd))  *  PWDOT 
BTERM2  «  Z1  *  V(I)  *((BNEW/Z1)**R2) 

C 

THETA2  =  BETAESTd)  *  DZ3SUB  /  Z3 
C 

BETADOTd)  =  AM2*BTERM1  -  A2*BTERM2 
C 

BETAESTd)  =  BETAd)  +  BETADOTd)*DTSUB  +  THETA2 
C 

ZD  =  ZD  +  BETAESTd)  *  Ud)  *  FACT 
CBIAXI  BSSUMl  =  BSSUMl  +  BETADOTd)  *  BETADOTd)  *  FACT 
CBIAXI  BSSUM2  =  BSSUM2  +  BETAESTd)  *  BETADOTd)  *  FACT 

44  CONTINUE 
C 

C  SPECIAL  TRAP  FOR  BETA  SNAP  REVERSALS  AT  HIGH  STRAIN  RATES 
C 

ZTOT  =  ZD  +  ZIEST 
IF(ZTOT.LT.O)THLT4 
BSCALE  *  ZIEST  /  ABS(ZD) 

ZD  =  ZD*BSCALE*0^ 

DO  45 1  =  1^ 

45  BETAESTd)*  BETAESTd) ‘BSCALE  *0.5 
ENDIF 

C 

BERROR  =  DABS(BNEW-BOLD2) 

BMAX  =  DMAXl(BNEW,BOLD2) 

IF(  BMAX  .GT.  TOLER) 
tc  BERROR  =  BERROR  /  BMAX 
C 

BOLD2  =  BNEW 
C 

IF  (BERROR  GT.  TOLER)  THEN 
IBCOUNT  =  IBCOUNT  +  1 
IF  (IBCOUNT  .GT.  5)  THEN 
ISUB  =  ISUB  *  2 
IF(  ISUB  .GT.  64)  IDBUG  =  1 
IF(  ISUB  .GT.  128)  THEN 
WRITE  (6,‘)  ■  BETA  REFUSES  TO  CONVERGE  ’ 

STOP 

ELSE 

GOTO  100 
ENDIF 
ENDIF 
CX)T041 
ENDIF 
C 

CBIAXI  IF(BSSUM1  .LT.  l.OE-12)  BSSUMl  =  l.E-12 
CBIAXI  VVDOT  =  BSSUM2  /  BSSUMl 
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CBIAXI  IF(ABS(VVDOT).LE.l.)  THEN 
CBIAXI  THETA  =  ACC6(VVDOT) 

CBIAXI  ELSEIF(  WDOT  .GT.  1 .)  THEN 

CBIAXI  THETA  =  0. 

CBIAXI  ELSE 

CBIAXI  THETA  =  3.141592654 

CBIAXI  ENDIF 
C 

CBIAXI  ALDOT  =  AM2  *  (ALPHA1*SIN(THETA)-ALPHA)  ‘  PWDOT 
C 

CBIAXI  ALEST  =  ALPHA  +  ALDOT  ♦  DTSUB 
C 

ZTERMl  =  Zl*(ABS(ZIEST-Z2)/Zir*Rl 
C 

THETAl  =  DZ2SUB 
C 

CBIAXI  ZIDOT  =  AM1*{Z1 + ALEST*Z3-ZIEST)*PWDOT 
CBIAXI  &  -  A1*ZTERM1 

C 

ZIDOT  =  AMI '(Zl-ZIESTIT  WDOT 
&  -  Al'ZTERMl 

C 

C  COMPUTE  UPDATED  ZTOT  VALUE 
C 

ZIEST  =  ZI  +  ZIDOT*DTSUB  +  THETAl 
C 

C  NOW  INVESTIGATE  CONVERGENCE 
C 

SERROR  =  DABSISIGEST-SIGNEW) 

SIGMAX  =  DMAXKSIGNEW^IGEST) 

1F(SIGMAX.GT.  TOLER) 

&  SERROR  =  SERROR  /  SIGMAX 
C 

DEPERR  =  DABS(DEIEFF-DOLD) 

DEPMAX  =  DMAXl{DElEFF,DOLD) 

IF{  DEPMAX  .GT.  TOLER) 

&  DEPERR  =  DEPERR  /  DEPMAX 
C 

ZERROR  =  DABS(ZIEST-ZOLD) 

ZMAX  =  DMAXl(ZIEST,ZOLD) 

IF(  ZMAX  .GT.  TOLER) 

&  ZERROR  =  ZERROR  /  ZMAX 
C 

ERROR  =  SERROR  +  DEPERR  +  BERROR  +  ZERROR 
C 

IFdDBUG.EQ.DTHEN 

WRITE(6,*) '  TIME  KTR  ICOUNT  AND  ISUB  ’JIME.  KTR,ICOUNT,ISUB 
WRITE(6,*) '  STRESS ' .  STRESS 

WR1TE(6,*)  'SIG  EST  NEW  SERR  SIGEST,SIGNEW,SERROR 
WRITE(6,*)  'DEP  OLD  NEW  EPERR DOLD,DEIEFF,DEPERR 
WR1TE(6,*)  'B  OLD  NEW  BERR  BOLD,BNEW,BERROR 
WRITE(6,*)  Z  OLD  NEW  ZERR ZOLD^IEST,ZERROR 
WRITE(6,*) 


73 


uuu  uuu  uuu  u  uuu  u  uuu  uuu  u  u  u  uuu 


ENDIF 


CHECK  FOR  CONVERGENCE 

IF(  ERROR.GT.  TOLER)  THEN 

NONCONVERGENT  SOLUTION  ARRIVES  HERE 

SIGEST  =  SIGOLD*(l. -RELAX)  +  SIGNEW  *  (RELAX) 
SIGOLD  =  SIGEST 
DOLD  =  DEIEFF 
BOLD  =  BNEW 
ZOLP  =  ZIEST 

UPDATE  AND  CHECK  CONVERGENCE  COUNT 
ICOUNT=ICOUNT+l 
IF  (ICOUT'JT  .GT.30)  THEN 
MAKE  ADDITIONAL  TIME-STEP  CUTS 
ISUB  =  ISUB  *  2 
IF  (ISUB  .EQ.  128)  IDBUG  =  1 
TERMINATE  NONCONVERGE  SOLUTIONS 
IF  (ISUB  .GT.  128)  THEN 

WRITE(6/)  ■  SOLUTION  REFUSED  TO  CONVERGE  ’ 
STOP 
ENDIF 

RESTART  ITERATION  SCHEME  WITH  NEW  TIME-STEP  CUT 
GOTO  100 
ENDIF 
GOTO  300 
ELSE 

CONVERGENCE  SOLUTION  ARRIVES  HERE 

ALPHA  =  ALEST 
ZI  =  ZIEST 
C 

DO  501  =  1,6 
EIN(I)  =EINEST(I) 

BETA(I)  =BETAEST(I) 

50  CONTINUE 
C 

ENDIF 
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C  END  OF  SUBTIME  LCX)PS 
C 

200  CONTINUE 
C 

SIGEFF  =  SlGEST 
BETAEFF  =  BNEW 
C 
C 

C  DEBUG  OUTPUT  STATEMENTS 
C 

CJOE2  IF  (NEL  .EQ.  1  .AND.  IPT  .EQ.  1  .AND.  KTR  .EQ.  INTER)  THEN 
C 

CJOE2  WRrrE(6  *)  ’  TIME  KTR  ICOUNT ISUB  IPT  ',TIME,  KTR,lCOUNT,ISUB,IPT 
CJOE2  WRITE(6,*) 

CJOE2  WRITE(6  •)  'SIG  EST  NEW  SERROR SIGEST,SIGNEW^ERROR 
CIOE2  WRITE<6,*)  'DEP  OLD  NEW  EPERROR  DOLD,DEIEFF,DEPERR 
CJOE2  WRITE<6,*)  'B  OLD  NEW  BERROR  BOLD,BNEW,BERROR 
CJOE2  WRITE(6,*)  ‘Z  OLD  NEW  2ERROR  *,  Z0LD^IEST;ZERR0R 
C]OE2  WRITE(6,*) 

CJOE2  WRITE(6,*) '  FOR  ELEMENT  AND  G  POINT  NEL,  IPT 
C|OE2  WRI're(6.*) 

CJOE2  WRITE(6,*) '  EINEST  EINEST 

CJOE2  WRITE{6,*)  ’  STRESS STRESS 

CJOE2  WRITE(6,*) '  DEVSIG  DEVSIG 

CJOE2  WRITE(6,*) '  STRAIN EPS 

CJOE2  WRITE(6,*)  ‘  BETAEST BETAEST 

CJOE2  WRITE(6  *) '  ZIEST,  DEIEFF,  ZD  ZIEST,  DEIEFF,  ZD 

C 

CJOE2  ENDIF 
C 

RETURN 

END 

C 

SUBROUTINE  STGET3(ARRAY) 

C 

IMPLICIT  DOUBLE  PRECISION  ( A-H,0-Z ) 

C 

REAL’S  ARRAY(60) 

C 

COMMON  /BPSTAT3/  EIN(6),  ZI,  SIGEFF,  BETAEFF,  EPEFF,  ALPHA, 

&  ZD,  BETA(6),  EPS0(6) 

C 

EIN(l)  =  ARRAY(I) 

EIN(2)  =  ARRAY(2) 

E1N(3)  =  ARRAY(3) 

EIN(4)  =  ARRAY(4) 

EIN(5)  =  ARRAY(5) 

EIN(6)  =  ARRAY(6) 

BETA(l)  =  ARRAY!?) 

BETA(2)  =  ARRAY(8) 

BETA(3)  =  ARRAY(9) 

BETA(4)  =  ARRAY(IO) 

BETA(5)  =  ARRAY(ll) 
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BETA(6)  =  ARRAY(12) 

ZI  =  ARRAY(13) 

SIGEFF  =  ARRAY(]4) 

EPEFF  rrARRAYdS) 

BETAEFF  =  ARRAYS  6) 

ALPHA  =ARRAY(17) 

ZD  =  ARRAYCtS) 

EPS0<1)  =  ARRAY(19) 

EPS0{2)  =  ARRAY(20) 

EPS0{3)  =  ARRAY(21) 

EPS0(4) »  ARRAY(22) 

FPS0(5)  =  ARRAY(23) 

EPS0(6)  =  ARRAY(24) 

C 

RETURN 

END 

C 

SUBROLTi  INE  STPUT3(ARRAY) 

C 

IMPLICIT  DOUBLE  PRECISION  (  A-H,0-Z  ) 

C 

REAL'S  ARRAY(60) 

C 

COMMON  /BPSTAT3/  EIN(6),  ZI,  SIGEFF,  BETAEFF.  EPEFF,  ALPHA, 
&  ZD,  BETA(6),  EPS0(6) 

C 

ARRAY(l)  =EIN(1) 

ARRAY(2)  =  EIN{2) 

ARRAY(3)  =EIN(3) 

ARRAY(4)  =  EIN(4) 

ARRAY(5)  =EIN(5) 

ARRAY(6)  =  EIN(5) 

ARRAY(7)  =BETA(1) 

ARRAYIS)  =  BETA(2) 

ARRAY(9)  =BETA(3) 

ARRAY(IO)  =  BETA(4> 

ARRAY(n)=  BETA(5) 

ARRAY(12)  =  BETA(6) 

ARRAY(13)  =  ZI 
ARRAY(14)  =  SIGEFF 
ARRAY!  15)  =  EPEFF 
ARRAY(]6)  =  BETAEFF 
ARRAY(17)  =  ALPHA 
ARRAY(18)  =  ZD 
ARRAY(19)  =  EPS0(1) 

ARRAY(20)  =  EPS0(2) 

ARRAY(21)  =  EPS0{3) 

ARRAY(22)  =  EPS0(4) 

ARRAY(23)  =  EPS(K5) 

ARRAY{24)  =  EI>S0(6) 

C 

RETURN 

END 


C 
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Appendix  C 


ADINAIN3.0  Input  File  for  Verification  Test  Case  1 
(file  vcasel.in) 


FEPROGRAM  ADINA 
DATABASE  CREATE  nLEUNIT=l 

HEADING  Test  case  1  of  Directional  Bodner-Partom  with  Beta21s’ 

MASTER  IDOF=100111  REACTIONS=YES  MTarK=500  NSTEP=48  DT=1.0, 
MODEX=EX  IRINT=48 
ANALYSIS  TYPE=STATIC 

ITERATION  METHOD=FULL  NEWTON  LINE-SEARCH=YES 
TOLERANCE  TYPE=ED  DNORM=l.e-5 

*  Displacement  Loading  with  time  (strain  rate  =  833.3e-6/s) 
TIMEFUNCTION  1 

0.0  0.0 

48.0  0.04 

*  Temperature  with  time 
TlMEFUCnON  2 

0.0  25.0 
48.0  25.0 

PRINTOUT  VOLUME=M  AX 

PORTHOLE  VOLUME=MAX  PORMATTED^YES 

* 

SYSTEM  1  CARTESIAN 

COORDINATES 

ENTRIES  NODE  ZL  YL 

1  O.lOOE+1  O.lOOE+1 

2  O.lOOE+1  O.OOOE-0 

3  O.OOOE-0  O.OOOE-0 

4  O.OOOE-0  O.lOOE+1 


•  Bodner  Partom  Test  Material  for  b21S 

*  with  no  Thermal  Expansion 

MATERIAL  1  USER  CTI]=0.(X)  CT12=1600.  CT13=3.  Cri4=l.c4, 
CT15=0.34  CT16=5.8c5  CT17=L37c4  CT18=0.0  , 

SCP1»0.75  SCP2=0.0001  SCP3=1  SCP4=1  , 

XINTER=10  TREF=25. 

*Temp(C)Expans(1/C) ElastidMPa)  n  zOd/s)  z3{MPa) m2{l/MPa) 

23.000  O.OOOOe^  112000.  4.8000  1550.0  100.0  0  350 

260.00  O.OOOOc-6  108030.  3.5000  1300.0  300.0  0.350 

315.00  O.OOOOe-6  106130.  3.0540  1250.4  390.0  1.502 

365.00  O.OOOOc-6  104090.  2.6490  1205.4  500.0  2.549 

415.00  O.OOOOe-6  101740.  2.2430  1160.4  660.0  3,597 

465.00  O.OOOOe-6  99085.  1.8380  1115.3  960.0  4.644 
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482.00  O.OOOOe-6 
500.00  0.0000e^ 
525.00  O.OOOOe-6 
550.00  O.OOOOe-6 
575.00  O.OOOOe-6 
600.00  O.OOOOe-6 
650.00  O.OOOOe-6 
760.00  O.OOOOe-6 
815.00  O.OOOOe-6 
900.00  O.OOOOe-6 


98113.  1.7000  1100.0  1100.0  5.000 

97045.  1.5000  1089.3  1300.0  5.763 

95497.  1.2800  1074.4  1670.0  6.822 

93873.  1.1000  1059.5  2100.0  7.881 

92172.  0.9700  1044.6  2600.0  8.941 

90395.  0.8200  1029.8  3700.0  10.000 

86612.  0.7400  1000.0  3800.0  10.000 

77216.  0.5800  600.0  4000.0  15.000 

71964.  0.5500  300.0  4100.0  30.000 

63122.  0.5500  300.0  4300.0  30.000 


EGROUP 1  TWODSOLID  SUBTYPE=AX1SYMMETR1C  MATER1AL=1 
GSURFACE  Nl=l,N2=2,N3=3,N4=4,ELl=l,EL2=:l,NODES=4 

BOUNDARIES  IDOF=110111  TYPE=NODES 
2 

BOUNDARIES  IDOF=llllll  TYPE=NODES 

3 

BOUNDARIES  IDOF=101111  TYPE=NODES 

4 

INITIAL  TEMPERATURES 

1  25.  STEP1T04  25. 

• 

LOADS  TEMPERATURE 
1 1.2  0STEP1T041.2  0 

LOADS  DISPLACEMENTS 
131.01 

2  31.01 

* 

*  ADINA  INPUT  HLE  CREATE 

ADINA 

END 
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Appendix  D 


Mathematica  Results  for  Verification  Test  Case  1 
(file  vcasel.math) 


Strain 

Stress  (MPa) 

Strain 

Stress  (MPa) 

0.0000000 

0.0000 

0.031665 

1146.7 

0.0008333 

93.330 

0.032499 

1146.7 

0.0016666 

186.66 

0.033332 

1146.7 

0.0024999 

279.99 

0.0^4165 

1146.7 

0.0033332 

373.32 

0.034999 

1146.7 

0.0041665 

466.65 

0.035832 

1146.7 

0.0049998 

559.98 

0.036665 

1146.7 

0.0058331 

653.31 

0.037499 

1146.7 

0.0066664 

746.64 

0.038332 

1146.7 

0.0074997 

839.97 

0.039165 

1146.7 

0.0083330 

933.30 

0.039998 

1146.7 

0.0091663 

1026.6 

0.0099996 

1083.9 

0.010833 

1098.6 

0.011666 

1110.4 

0.012500 

1119.6 

0.013333 

1126.7 

0.014166 

1.132.0 

0.014999 

1136.0 

0.015833 

1138.9 

0.016666 

1141.1 

0.017499 

1142.7 

0.018333 

1143.8 

0.019166 

1144.6 

0.019999 

1145.2 

0.020832 

1145.6 

0.021666 

1145.9 

0.022499 

1146.2 

0.023332 

1146.3 

0.024166 

1146.4 

0.024999 

1146.5 

0.025832 

1146.6 

0.026666 

1146.6 

0.027499 

1146.7 

0.028332 

1146.7 

0.029166 

1146.7 

0.029999 

1146.7 

0.030832 

1146.7 
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Appendix  E 


Mathematica  Results  for  Verification  Test  Case  2 
(file  vcase2.math) 


Strain 

Stress  (MPa) 

Strain 

Stress  (MPa) 

0.0000000 

0.0000 

0.031665 

140.38 

0.0008333 

72.106 

0.032499 

140.38 

0.0016666 

121.44 

0.033332 

140.38 

0.0024999 

138.00 

0.034165 

140.38 

0.0033332 

140.14 

0.034999 

140.38 

0.0041665 

140.36 

0.035832 

140.38 

0.0049998 

140.38 

0.036665 

140.38 

0.0058331 

140.38 

0.037499 

140.38 

0.0066664 

140.38 

0.038332 

140.38 

0.0074997 

140.38 

0.039165 

140.38 

0.0083330 

140.38 

0.039998 

140.38 

0.0091663 

140.38 

0.0099996 

140.38 

0.010833 

140.38 

0.011666 

140.38 

0.012500 

140.38 

0.013333 

140.38 

0.014166 

140.38 

0.014999 

140.38 

0.015833 

140.38 

0.016666 

140.38 

0.017499 

140.38 

0.018333 

140.38 

0.019166 

140.38 

0.019999 

140.38 

0.020832 

140.38 

0.021666 

140.38 

0.022499 

140.38 

0.023332 

140.38 

0.024166 

140.38 

0.024999 

140.38 

0.025832 

140.38 

0.026666 

140.38 

0.027499 

140.38 

0.028332 

140.38 

0.029166 

140.38 

0.029999 

140.38 

0.030832 

140.38 
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Appendix  F 


Mathematica  Results  for  Verification  Test  Case  3 
(file  vcaseS.math) 


Strain 

Stress  (MPa) 

Strain 

Stress  (MPa) 

O.OOOOOOO 

0.0000 

0.031665 

140.38 

0.0008333 

72.106 

0.032499 

140.38 

0.0016666 

121.44 

0.033332 

140.38 

0.0024999 

138.00 

0.034165 

140.38 

0.0033332 

140.14 

0.034999 

140.38 

0.0041665 

140.36 

0.035832 

140.38 

0.0049998 

140.38 

0.036665 

140,38 

0.0058331 

140.38 

0.037499 

140.38 

0.0066664 

140.38 

0.038332 

140.38 

0.0074997 

140..38 

0.039165 

140.38 

0.0083330 

140.38 

0.039998 

140.38 

0.0091663 

140.38 

0.0099996 

140.38 

0.010833 

140.38 

0.011666 

140.38 

0.012500 

140.38 

0.013333 

140.38 

0.014166 

140.38 

0.014999 

140.38 

0.015833 

140.38 

0.016666 

140.38 

0.017499 

140.38 

0  018333 

140.38 

0.019166 

140.38 

0.019999 

140.38 

0.020832 

140.38 

0.021666 

140.38 

0.022499 

140.38 

0.023332 

140.38 

0.024166 

140.38 

0.024999 

140.38 

0.025832 

140.38 

0,026666 

140.38 

0.027499 

140.38 

0.028332 

140.38 

0.029166 

140.38 

0.029999 

140.38 

0,030832 

140.38 
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Appendix  G 

Mathematica  Results  for  Verification  Test  Case  4  (file  vcase4.math) 


Strain 

Stress 

(MPa) 

Temp  (°C) 

Strain 

Stress 

(MPa) 

Temp  (°C) 

0.00000000 

0.0000 

25.000 

0.032500 

787.03 

482.00 

0.00083333 

93.333 

25.000 

0.033333 

787.03 

482.00 

0.0016667 

186.67 

25.000 

0.034167 

787.03 

482.00 

0.0025000 

280.00 

25.000 

0.035000 

787.03 

482.00 

0.0033333 

373.33 

25.000 

0.035833 

787.03 

482.00 

0.0041667 

466.67 

25.000 

0.036667 

787.03 

482.00 

0.0050000 

560.00 

25.000 

0.037500 

787.03 

482.00 

0.0058333 

653.33 

25.000 

0.038333 

787.03 

482.00 

0.0066667 

746.67 

25.000 

0.039167 

787.03 

482.00 

0.0075000 

840.00 

25.000 

0.040000 

787.03 

482.00 

0.0083333 

933.33 

25.000 

0.040833 

787.03 

482.00 

0.0091667 

1026.7 

25.000 

0.041667 

787.03 

482.00 

0.0100000 

1083.9 

25.000 

0.042500 

787.03 

482.00 

0.010833 

1098.6 

25.000 

0.043333 

787.03 

482.00 

0.011667 

1110.4 

25.000 

0.044167 

787.03 

482.00 

0.012500 

1119.6 

25.000 

0.045000 

787.03 

482.00 

0.013333 

1126.7 

25.000 

0.045833 

787.03 

482.00 

0.014167 

1132.0 

25.000 

0.046667 

787.03 

482.00 

0.015000 

1136.0 

25.000 

0.047500 

787.03 

482.00 

0.015833 

1138.9 

25.000 

0.048333 

787.03 

482.00 

0.016667 

1141,1 

25.000 

0.049167 

787.03 

482.00 

0.017500 

1142.7 

25.000 

0.050000 

787.03 

482.00 

0.018333 

1143.8 

25.000 

0.050833 

821.57 

443.90 

0.019167 

1144.6 

25.000 

0.051667 

843.77 

405.80 

0.020000 

1145.2 

25.000 

0.052500 

879.05 

367.80 

0.020833 

1122.6 

63.100 

0.053333 

913.54 

329.70 

0.021667 

1097.2 

101.20 

0.054167 

945.70 

291.60 

0.022500 

1070.3 

139.30 

0.055000 

975.69 

253.50 

0.023333 

1041.9 

177.30 

0.055833 

1008.1 

215.40 

0.024167 

1011.7 

215.40 

0.056667 

1038.7 

177.30 

0.025000 

979.53 

253.50 

0.057500 

1067.4 

139.30 

0.025833 

949.64 

291.60 

0.058333 

1094.6 

101.20 

0.026667 

917.28 

329.70 

0.059167 

1120.3 

63.100 

0.027500 

882.53 

367.80 

0.060000 

1144.7 

25.000 

0.028333 

850.18 

405.80 

0.060833 

1146.7 

25.000 

0.029167 

825.88 

443.90 

0.061667 

1146.7 

25.000 

0.030000 

790.59 

482.00 

0.062500 

1146.7 

25.000 

0.030833 

787.05 

482.00 

0.063333 

1146.7 

25.000 

0.031667 

787.03 

482.00 

0.064167 

1146.7 

25.000 
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Appendix  H 

Mathematica  Results  for  Verification  Test  Case  5  (file  vcaseS.math) 


Strain 

Stress 

(MPa) 

Temp  (®C) 

Strain 

Stress 

(MPa) 

Temp  i°C) 

0.0000 

0.0000 

650.00 

0.016250 

130.64 

760.00 

0.00041667 

36.095 

650.00 

0.016667 

129.59 

760.00 

0.00083333 

72.190 

650.00 

0.017083 

129.05 

760.00 

0.0012500 

104.07 

650.00 

0.017500 

128.76 

760.00 

0.0016667 

133.08 

650.00 

0.017917 

128.62 

760.00 

0.0020833 

162.60 

650.00 

0.018333 

128.54 

760.00 

0.0025000 

192.14 

650.00 

0.018750 

128.50 

760.00 

0.0029167 

221.16 

^50.00 

0.019167 

128.48 

760.00 

0.0033333 

249.09 

650.00 

0.019583 

128.47 

760.00 

0.0037500 

275.30 

650.00 

0.020000 

128.46 

760.00 

0.0041667 

299.15 

650.00 

0.020417 

128.46 

760.00 

0.0045833 

320.03 

650.00 

0.020833 

128.46 

760.00 

0.0050000 

337.52 

650.00 

0.021250 

128.46 

760.00 

0.0054167 

351.44 

650.00 

0.021667 

128.46 

760.00 

0.0058333 

361.94 

650.00 

0.022083 

128.46 

760.00 

0.0062500 

369.48 

650.00 

0.022500 

128.46 

760.00 

0.0066667 

374.65 

650.00 

0.022917 

128.46 

760.00 

0.0070833 

378.08 

650.00 

0.023333 

128.46 

760.00 

0.0075000 

380.30 

650.00 

0.023750 

128.46 

760.00 

0.0079167 

381.71 

650.00 

0.024167 

128.46 

760.00 

0.0083333 

382.60 

650.00 

0.024583 

128.46 

760.00 

0.0087500 

383.15 

650.00 

0.025000 

128.46 

760.00 

0.0091667 

383.49 

650.00 

0.025417 

134.59 

750.80 

0.0095833 

383.70 

650.00 

0.025833 

144.80 

741.70 

0.0100000 

383.83 

650.00 

0.026250 

157.42 

732.50 

0.010417 

374.19 

659.20 

0.026667 

171.74 

723.30 

0.010833 

357.33 

668.30 

0.027083 

187.43 

714.20 

0.011250 

336.61 

677.50 

0.027500 

204.29 

705.00 

0,011667 

314.21 

686.70 

0.027917 

222.20 

695.80 

0.012083 

291.34 

695.80 

0.028333 

241.07 

686.70 

0.012500 

268.60 

705.00 

0.028750 

260.85 

677.50 

0.012917 

246.28 

714.20 

0.029167 

281.50 

668.30 

0.013333 

224.60 

723.30 

0.029583 

302.96 

659.20 

0.013750 

203.67 

732.50 

0.030000 

325.21 

650.00 

0.014167 

183.61 

741.70 

0.030417 

341.88 

650.00 

0.014583 

164.50 

750.80 

0.030833 

354.79 

650.00 

0.015000 

146.43 

760.00 

0.031250 

364.38 

650.00 

0.015417 

136.98 

760.00 

0.031667 

371.18 

650.00 

0.015833 

132.73 

760.00 

0.032083 

375.79 

650.00 
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Strain 

Stress 

(MPa) 

Temp  (°C) 

0.032500 

378.83 

650.00 

0.032917 

380.78 

650.00 

0.033333 

382.01 

650.00 

0.033750 

382.78 

650.00 

0.034167 

383.26 

650.00 

0.034583 

383.56 

650.00 

0.035000 

383.74 

650.00 

0.035417 

383.86 

650.00 

0.035833 

383.93 

650.00 

0.036250 

383.97 

650.00 

0.036667 

383.99 

650.00 

0.037083 

384.01 

650.00 

0.037500 

384.02 

650.00 

0.037917 

384.03 

650.00 

0.038333 

384.03 

650.00 

0.038750 

384.03 

650.00 

0.039167 

384.04 

650.00 

0.039583 

384.04 

650.00 

0.040000 

384.04 

650.00 

4 

t 
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